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ABSTRACT 


The  area  comprising  southwestern  New  Mexico  and 
southeastern  Arizona  has  experienced  a  complex  tectonic 
history.  In  particular,  the  period  of  time  from  the  late 
Cretaceous  to  the  present  has  brought  varying  degrees  of 
compression,  magmatic  activity,  uplift  and  extension.  Three 
major  provinces  that  developed  as  a  result  of  this  tectonism 
are  the  Colorado  Plateau,  Basin  and  Range  and  Rio  Grande 
Rift. 

The  Colorado  Plateau  is  a  region  which  is  uplifted  and 
relatively  undeformed  with  respect  to  surrounding  provinces. 
It  is  characterized  by  gently  dipping  strata  that  have 
undergone  minor  folding  and  warping,  volcanism,  and 
epeirogenic  uplift  during  Cenozoic  time.  To  the  west,  south 
and  southeast,  the  Colorado  Plateau  is  bounded  by  the  Basin 
and  Range  and  Rio  Grande  Rift  extensional  provinces.  They 
have  undergone  extensive  deformation  and  volcanic  activity 
during  the  past  40  Ma,  with  signs  of  active  tectonism 
continuing  to  the  present. 

The  transition  zones  between  the  Colorado  Plateau,  Rio 
Grande  Rift  and  Basin  and  Range  provinces  meet  in 
southeastern  Arizona  and  southwestern  New  Mexico.  There  is 
no  physiographic  boundary  nor  is  there  a  distinctive  change 
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in  structural  style  observed  in  this  region.  Instead,  the 
surface  is  covered  by  the  mid-Tertiary  Mogollon-Datil 
volcanic  field,  which  hides  underlying  structures  that  might 
aid  in  determining  the  nature  of  the  transition.  Therefore, 
the  use  of  geophysical  techniques  in  conjunction  with  known 
geology  must  be  employed  to  determine  the  present-day  state 
of  the  lithosphere,  as  well  as  its  tectonic  history. 

A  combined  analysis  of  seismic  refraction  and  gravity 
data  has  identified  a  major  upper  crustal  low-velocity,  low- 
density  body  in  this  region.  This  body  is  associated  with  a 
thickening  of  the  crust  where  it  is  present.  It  underlies 
the  Mogollon-Datil  volcanic  field,  but  extends  to  the  north 
through  the  Plains  of  St.  Augustine  and  under  the  Datil 
Mountains.  Because  of  its  proximity  to  the  volcanic  field, 
this  body  is  interpreted  as  a  massive  plutonic  complex 
nearly  170  km  long  by  130  km  wide  in  the  study  area.  It 
presumably  is  the  source  for  the  volcanics,  implying  that 
this  complex  is  mid-Tertiary  in  age.  When  correlated  with 
long  wavelength  gravity  trends,  the  complex  extends 
northwest  into  Arizona,  and  is  approximately  230  km  by  155 
km  in  extent. 

Beneath  this  plutonic  complex,  seismic  and  gravity  data 
indicate  that  the  Moho  deepens  markedly  from  32  to  37  km 
about  60  km  north  of  Silver  City.  This  is  identified  as  the 
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present-day  structural  boundary  between  the  Colorado  Plateau 
and  the  Basin  and  Range  and  Rio  Grande  Rift  provinces. 
Compressional  wave  arrivals  and  gravity  analysis  demonstrate 
anomalously  slow  velocities  and  low  density  in  the  upper 
mantle  along  both  profiles  examined  in  this  study.  This 
supports  heat  flow  studies  which  conclude  that  the  upper 
mantle  and  lower  crust  contain  a  thermal  anomaly.  Although 
high  heat  flow  values  are  observed  in  both  the  Basin  and 
Range  and  Rio  Grande  Rift,  the  highest  heat  flow  values  are 
spatially  associated  with  the  transition  zone.  Because  of 
the  extensive  volcanism  and  evidence  for  geothermal  activity 
throughout  the  Cenozoic,  thf  thermal  anomaly  is  interpreted 
to  be  present  over  the  same  time  interval.  A  two- 
dimensional  finite  difference  model  has  therefore  been 
derived  and  a  computer  code  generated  which  maps  the 
progression  of  a  temperature  pulse  over  time.  Surface  heat 
flow  and  isostatic  uplift  are  also  calculated  by  this 
program.  These  values  are  compared  to  observed  values  of 
heat  flow  and  elevation  in  Arizona  and  New  Mexico  to 
constrain  both  the  magnitude  and  width  of  the  thermal  pulse. 

Results  of  the  analysis  are  consistent  with  thermal 
anomalies  present  at  both  the  base  of  the  crust  and  the  base 
of  the  lithosphere.  The  deep  anomaly  produces  a  steady- 
state  solution  after  nearly  28  Ma,  producing  the  bulk  of 
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isostatic  uplift  but  making  little  contribution  to  surface 
heat  flow.  The  shallow  anomaly  becomes  steady-state  after 
5.4  Ma,  and  is  responsible  for  present-day  surface  heat  flow 
readings.  It  contributes  a  minor  component  of  uplift  to  the 
region.  The  time  span  over  which  these  thermal 
perturbations  operate  implies  that  the  deeper  anomaly  is  the 
cause  for  much  of  the  volcanic  activity  since  Laramide  time, 
and  that  the  shallow  anomaly  is  associated  with  recent 
extension  and  basaltic  volcanism. 
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INTRODUCTION 


The  Colorado  Plateau  is  one  of  the  major  tectonic  and 
physiographic  provinces  of  western  North  America.  It  is 
approximately  3.5  x  10^  km^  in  area,  and  is  centered  near 
the  common  boundaries  of  New  Mexico,  Colorado,  Utah  and 
Arizona  (Figure  1) .  It  is  bounded  on  three  sides  by  the 
extensional  tectonic  provinces  of  the  Basin  and  Range  and 
Rio  Grande  Rift,  and  to  the  north  and  northeast  by  the 
Middle  and  Southern  Rockies,  respectively.  The  origin  of 
the  plateau,  its  late  Mesozoic  and  Cenozoic  evolution,  its 
deep  structure  and  its  relationship  to  adjacent  provinces 
pose  questions  of  broad  interest  to  the  earth  science 
community. 

The  interior  of  the  Colorado  Plateau  appears  to  be 
relatively  unaffected  by  the  extension  and  associated 
volcanism  that  border  it  on  three  sides.  Deformation  within 
the  plateau  during  and  since  Laramide  tectonism  has  been 
restricted  to  gentle  warping  and  minor  faulting.  Cenozoic 
volcanic  activity  has  been  less  extensive  on  the  plateau 
than  in  most  other  areas  of  the  western  United  States 
(Christiansen  and  Lipman,  1972;  Lipman  and  others,  1972), 
but  its  margins  have  been  loci  of  igneous  activity  (Morgan 
and  Swanberg,  1985) . 


Figure  l.  Map  of  the  western  United  States  showing  locations 
and  boundaries  of  major  tectonic  provinces. 
Modified  after  Thompson  and  Zoback  (1979) . 
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Despite  its  relative  tectonic  stability  and  lack  of 
magmatism,  the  Colorado  Plateau  has  been  epeirogenically 
uplifted  approximately  2  km  since  the  latest  Mesozoic,  based 
on  the  present  elevation  of  the  Cretaceous  (marine)  Mancos 
shale  (Morgan  and  Swanberg,  1985) .  Although  the  timing  of 
uplift  is  still  debated,  the  basic  geologic  history  of  the 
plateau  is  reasonably  well  documented  (Hunt,  1956) . 

The  margins  where  the  Colorado  Plateau  is  bounded  by 
the  Basin  and  Range  and  Rio  Grande  Rift  provinces  have  been 
tectonically  active  and  are  of  much  current  interest  (e.g. 
Keller  and  others,  1975;  Aldrich  and  Laughlin,  1984; 
Brumbaugh,  1987) .  These  margins  are  divided  into  three 
distinct  sectors:  western  (Utah  and  northwestern  Arizona) , 
southern  (central  Arizona)  and  southeastern  (southeast 
Arizona  and  southwest  New  Mexico) .  These  boundaries  have 
been  classically  drawn  on  a  physiographic  basis,  although  no 
physiographic  expression  exists  along  the  southeastern 
margin. 

The  basis  for  defining  province  boundaries  by 
physiographic  means  has  left  room  for  much  debate.  There 
exists  a  considerable  body  of  geophysical  information 
indicating  that  the  margins  are  transitional  (Keller  and 
others,  1975;  Thompson  and  Zoback,  1979,  Brumbaugh,  1987). 
Furthermore,  in  those  areas  where  the  physiographic 
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boundaries  are  not  well  defined,  the  placement  of  the 
transition  zone  is  still  in  question  (Brumbaugh,  1987) .  of 
particular  interest  to  this  study  is  the  transition  between 
the  Colorado  Plateau,  Basin  and  Range,  and  Rio  Grande  Rift 
provinces  in  southeastern  Arizona  and  southwestern  New 
Mexico  (Figure  2) . 

The  Mogollon-Datil  volcanic  field  (MDVF)  lies  within 
the  Colorado  Plateau-Basin  and  Range-Rio  Grande  Rift 
transition  zone  in  southwest  New  Mexico.  It  hides 
underlying  structures  that  may  aid  in  determining  the  extent 
and  nature  of  this  transition  (Elston  and  others,  1976a; 

Gish  and  others,  1981) .  A  careful  integrated  geophysical 
and  geological  analysis  is  required  to  determine  the 
validity  of  previous  models  and  determine  the  deep  structure 
of  this  region. 

A  gravity  study  by  Daggett  and  others  (1986)  delineated 
two  regional  gravity  anomalies  in  southwestern  New  Mexico. 

A  positive  anomaly  occurs  along  the  Rio  Grande  Rift  axis 
near  El  Paso,  and  a  negative  anomaly  corresponds  to  the 
MDVF.  Their  crustal  model  indicates  crustal  thickening 
under  the  MDVF.  Sinno  and  others  (1986)  also  interpret  an 
increase  in  crustal  thickness  from  the  Rio  Grande  Rift  into 
the  MDVF  based  on  seismic  data.  Gish  and  others  (1981) 
concluded  that  crustal  thickening  occurs  in  the  MDVF  with 


Figure  2.  Location  of  study  area. 
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respect  to  southern  Arizona,  primarily  caused  by  a  five  km 
thick  low-density  surface  layer  composed  of  intermediate  to 
felsic  volcanic  rocks.  Further  north,  however,  Jaksha 
(1982)  suggests  that  little  crustal  thickening  occurs  across 
the  field  between  the  Rio  Grande  and  Morenci,  Arizona. 

Elston  (1976),  Rhodes  (1976),  Elston  and  others 
(1976a),  Krohn  (1976)  and  Coney  (1976a)  agree  that  a  shallow 
Plutonic  complex  underlies  the  high-silica  rhyolite 
extrusives  of  the  Mogollon  plateau.  Gravity  anomalies  were 
used  to  infer  that  this  complex  is  presumably  less  dense 
than  the  material  it  displaced  (Coney,  1976a;  Krohn,  1976), 
but  the  boundaries  and  thickness  of  the  body  were  not 
determined. 

The  extrusives  of  the  MDVF  have  completely  buried  the 
material  which  previously  formed  the  surface.  To  the  south 
of  the  volcanics.  Basin  and  Range  structures  predominate. 

The  Rio  Grande  Rift  lies  immediately  to  the  east,  and 
Colorado  Plateau  structural  features  predominate  to  the 
northwest.  However,  no  physiographic  boundary  between  the 
provinces  is  discernible  within  the  MDVF.  Thus,  a  major 
objective  of  this  study  is  to  use  a  combination  of 
geophysical  techniques  and  regional  geologic  data  to  study 
the  nature  and  history  of  the  transition  zone  between  the 
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southeastern  Colorado  Plateau  and  the  Rio  Grande  Rift  and 
Basin  and  Range  province. 

The  present  day  state  of  the  lithosphere  is  a 
culmination  of  preceding  tectonic  events.  Understanding  the 
thermal  history  of  the  region  is  therefore  a  second  goal  of 
this  study.  Post  mid-Tertiary  history  in  southwest  New 
Mexico  and  southeast  Arizona  is  dominated  by  block  faulting 
and  basaltic  volcanism  associated  with  Rio  Grande  Rift  and 
Basin  and  Range  extension.  The  volcanism  within  the  rift, 
however,  has  been  of  relatively  minor  volume  when  compared 
to  the  similarly-dimensioned  East  African  Rift  in  Kenya 
(Baldridge  and  others,  1984).  The  presence  of  hot  springs 
and  late  Cenozoic  volcanism  along  with  observed  geophysical 
data  indicate  that  elevated  temperatures  are  present  in  the 
lower  crust  and  upper  mantle  (Decker  and  Smithson,  1975; 
Seager  and  Morgan,  1979;  Gish  and  others,  1981;  Morgan  and 
others,  1986;  Sinno  and  others,  1986).  Heat  flow  data  are 
elevated  near  the  late  Cenozoic  volcanic  centers  and  near 
deep  basin-bounding  faults  in  the  region,  but  are  of  more 
intermediate  values  elsewhere.  However,  the  heat  flow  data 
are  complicated  in  the  area  by  extensive  groundwater 
circulation,  and  therefore  must  be  interpreted  with  caution 
(Morgan  and  others,  1986) . 
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Within  this  context,  a  broad  zone  of  elevated  heat  flow 
values  is  observed  near  the  boundary  of  the  Colorado  Plateau 
in  southeastern  Arizona  and  south  and  central  New  Mexico 
(Swanberg,  1979) .  Combined  with  the  crustal  uplift  in  the 
region  (Morgan  and  Swanberg,  1985) ,  which  is  at  least 
partially  due  to  thermal  expansion,  a  long-standing  thermal 
anomaly  is  inferred  at  the  base  of  the  lithosphere.  A  major 
goal  of  this  study,  therefore,  is  the  development  of  a  two- 
dimensional  finite  difference  model  in  which  the  evolution 
of  temperature,  heat  flow  and  isostatic  uplift  is  calculated 
over  time.  The  results  of  this  program  are  compared  with 
observed  heat  flow  and  elevation  for  the  purpose  of  defining 
the  thermal  history  of  the  region. 


GEOLOGIC  HISTORY 


Precambrian 

Precambrian  basement  is  exposed  in  scattered  locations 
throughout  the  region  (Dane  and  Bachman,  1965;  Wilson  and 
others,  1969) .  It  is  composed  primarily  of  crystalline 
schists,  gneisses,  and  granites  to  the  west;  and  slightly 
younger  Precambrian  sediments  in  the  eastern  part  of  the 
region  (Ross  and  Ross,  1986;  Seager  and  Mack,  1986).  The 
basement  is  usually  unconformably  overlain  by  Paleozoic 
rocks . 

The  Precambrian  basement  in  southern  Arizona  and  New 
Mexico  consists  of  a  complex  assemblage  of  rocks  ranging  in 
age  from  1.68  to  1.8  Ga  using  U-Pb  zircon  dates  (Figure  3; 
Condie,  1986).  These  ages  are  supported  by  Nd  isotope 
studies  (Bennett  and  DePaolo,  1987) .  The  rocks  show 
evidence  of  tectonic  activity  prior  to  and  during 
emplacement  (Grambling  and  others,  1988;  Karlstrom  and 
Bowring,  1988) .  The  structural  grain  in  the  Arizona 
transition  zone  strikes  northeast,  becoming  more  easterly  in 
central  New  Mexico.  Anorogenic  granites  (1.4  to  1.5  Ga) 
intrude  all  rocks  in  the  southwestern  United  States,  and  are 
especially  prevalent  in  southeast  Arizona  and  southern  New 
Mexico  (Condie,  1982) .  The  region  experienced  an  anorogenic 
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Figure  3 .  Precambrian  provinces  of  North  America .  After 
Condie  (1986) . 
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volcano-plutonic  event  of  middle  Proterozoic  age  (1.1  Ga) , 
as  recorded  in  the  Franklin  mountains  in  El  Paso  (Shannon, 
1989)  . 

Karlstrom  and  Bowring  (1988)  have  mapped  at  least  eight 
tectonostratigraphic  terranes  along  the  transition  zone  in 
Arizona.  The  easternmost  is  the  Pinal  block  adjacent  to 
Cenozoic  volcanic  cover  near  the  New  Mexico  border.  To  the 
east  of  Cenozoic  volcanic  cover,  six  tectonostratigraphic 
terranes  are  reported  in  New  Mexico  (Grambling  and  others, 
1988)  .  The  structural  grain  resulting  from  emplacement  of 
these  blocks  may  have  exerted  considerable  control  on  post- 
Precambrian  regional  tectonic  style. 

Paleozoic 

The  Paleozoic  section  in  southeastern  Arizona  reflects 
deposition  in  a  stable,  cratonal  shelf  environment  with  no 
highlands  nearby  (Coney,  1978a) .  Despite  this  indication  of 
regional  stability,  there  were  two  Paleozoic  erogenic  events 
in  southern  and  western  North  America:  1)  within  800  km  to 
the  west,  the  Cordilleran  miogeocline  experienced  the 
effects  of  the  Antler  and  Sonora  erogenic  events;  and  2)  the 
Appalachian-Ouachita-Marathon  orogenies,  resulting  from 
interactions  between  North  America  and  Africa-South  America, 
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came  within  500  km  to  the  east.  It  is  possible  that  these 
two  belts  may  have  intersected  somewhere  to  the  south  of 
Arizona,  possibly  within  1000  km.  Southeastern  Arizona  can 
therefore  be  thought  of  as  part  of  a  southwest-trending 
peninsular  extension  of  the  Paleozoic  North  American  craton. 

The  closure  of  the  lapetus  ocean  and  ancestral  Gulf  of 
Mexico  initiated  two  important  events  in  the  present  study 
area:  1)  the  Pennsylvanian  Ouachita-Marathon  orogeny,  which 

created  the  Ancestral  Rocky  Mountains  and  produced  faulting 
and  zones  of  weakened  crust  in  the  region  (Coney,  1978a, 
1978b;  Kluth,  1986;  Ross  and  Ross,  1986).  2)  although 

southeastern  Arizona  was  relatively  stable  throughout  this 
time,  southwestern  New  Mexico  lay  at  the  southwest  boundary 
of  a  large  basin-uplift  regime  that  extended  onto  the 
present  Colorado  Plateau  (Figure  4).  By  the  late  Early 
Permian,  the  entire  region  was  uplifted,  and  sedimentation 
was  mostly  nonmarine  to  the  end  of  the  period  (Ross  and 
Ross,  1986) . 


Mesozoic 

South  of  the  Colorado  Plateau,  little  information 
exists  from  the  early  Mesozoic  (Coney,  1978a;  Jahns  and 
others,  1978;  Seager  and  Mack,  1986).  However,  a  Jurassic 
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Figure  4.  Map  showing  Paleozoic  tectonic  features  of  the 

southwestern  United  States  and  northern  Chihuahua, 
Mexico.  After  Ross  and  Ross  (1986). 
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magmatic  zone  existed  in  Arizona  (Coney,  1978a) .  It  had  a 
northwest-southeast  trend,  parallel  to  and  just  south  of  the 
present  day  transition  zone  between  the  Colorado  Plateau  and 
Basin  and  Range  provinces.  It  continued  into  Mexico, 
passing  just  to  the  southwest  of  New  Mexico. 

This  zone  was  summarized  by  Hamilton  (1969) ,  Burchfiel 
and  Davis  (1972) ,  and  Coney  (1978a)  as  being  representative 
of  an  Andean-type  arc  setting,  implying  a  compressional 
regime  in  late  Triassic  to  early  Jurassic  time.  However, 
citing  modern  analogs  in  Cental  America  and  Kamchatka, 
Busby-Spera  (1988)  suggests  an  extensional,  graben-related 
volcanic  episode.  Sediments  from  volcanic  complexes  were 
not  shed  to  the  north  and  northeast,  due  to  the  presence  of 
a  low-lying  Precambrian-cored  ancestral  Mogollon  Highland, 
near  the  present  transition  zone  (Coney,  1978a) . 

The  Gulf  of  Mexico  began  to  open  during  Jurassic  time 
(Coney,  1978a) .  This  produced  two  events  of  note  in 
southern  Arizona  and  New  Mexico,  one  structural  and  one 
stratigraphic:  1)  as  the  gulf  opened,  a  major  transform 

fault  or  "megashear"  (Silver  and  Anderson,  1974)  may  have 
developed  on  the  southwest  edge  of  the  magmatic  arc, 
transporting  much  of  northern  Mexico  to  its  present 
location,  and  2)  the  opening  of  the  gulf  initiated  a  marine 
transgression  that  grew  out  of  the  gulf  into  western  North 
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America  (Coney,  1978a,  1978b) .  The  earliest  known  Mesozoic 
sediments  in  southern  New  Mexico  are  marine  Jurassic  rocks 
found  in  a  single  oil  test  well  near  Las  Cruces  (Thompson 
and  Bieberman,  1975) . 

Tectonism  and  volcanism  diminished  in  southern  Arizona 
and  New  Mexico  during  the  latter  part  of  the  Jurassic 
(Coney,  1978a,  1978b) .  During  the  Cretaceous,  the  Sevier 
orogeny  commenced  in  western  North  America  (Armstrong,  1968; 
Coney,  1978b) .  In  the  southwestern  United  States,  the 
volcanic  arc  was  reinstated  to  the  southwest  of  the  Jurassic 
arc  (Figure  3,  Coney,  1978a).  A  remnant  of  the  Mogollon 
Highland  persisted  as  the  Burro-Deming  Axis  from  the  West 
Potrillo  Mountains  to  the  northwest  into  Arizona  (Seager  and 
Mack,  1986) ,  but  the  rest  of  the  region  was  covered  by  a 
shallow  sea.  Near  the  arc,  significant  deformation  occurred 
in  southwest  Arizona  and  northwest  Sonora.  This  created  a 
highland  that  shed  sediments  towards  southwest  New  Mexico 
(Coney,  1978a) .  It  is  unclear  whether  this  deformation  was 
a  precursor  to  the  Laramide  compressional  event,  yet  it 
apparently  predates  Laramide  plutonism  in  the  region. 
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Laramide 

The  Laramide  orogeny  began  in  late  Cretaceous  time,  and 
involved  the  entire  crust  in  the  southwestern  United  States 
(Coney,  1978a,  1978b) .  In  western  North  America,  the 
volcanic  arc  broke  up  into  scattered  volcano-plutonic 
complexes.  Basement  cored  Laramide  uplifts  developed  where 
gaps  in  volcanism  occurred. 

In  southern  Arizona  and  New  Mexico,  Epis  and  Chapin 
(1975)  report  that  the  Laramide  event  began  with  widespread 
erosion  to  a  surface  of  little  relief.  This  was  followed  by 
uplift,  volcanism,  compression,  then  further  volcanism  and 
plutonism.  Coney  (1978a)  states  that  volcanism  progressed 
from  west  to  east  (into  Trans-Pecos  Texas),  then  swept 
westward  again.  This  ’’retrograde  motion"  of  volcanism  swept 
through  heated  crustal  material. 

Compressive  deformation  in  southern  Arizona  and  New 
Mexico  is  typical  of  that  found  throughout  the  western 
United  States  during  Laramide  time,  but  individual 
structures  are  smaller  (Seager  and  Mack,  1986) .  In  southern 
New  Mexico,  Laramide  compressive  features  are  displayed  as 
west-northwest  trending  asymmetric  block  uplifts  with 
complementary  basins.  The  asymmetry  is  found  in  the  style 
of  uplift,  where  the  northeast  boundary  generally  consists 
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of  reverse  and  thrust  fault  zones  that  dip  to  the  southwest 
under  the  uplift  (Figure  5) .  Laramide  basin  sediments  lap 
onto  the  unfaulted  southwest  edge  of  the  uplifts.  These 
structures  are  spatially  and  temporally  related  to 
structures  in  the  Rocky  Mountains  (Berg,  1962;  Tweto,  1975; 
Gather  and  Johnson,  1986) ,  thin-skinned  deformation  in  the 
Sierra  Madre  Oriental  (Coney,  1976b) ,  and  to  the  Laramide 
monoclines  of  the  Colorado  Plateau  (Kelley,  1955;  Coney, 
1976b;  Davis,  1978). 

The  thrust  faults  often  contain  enough  strike-slip 
motion  to  be  termed  oblique  faults  (Seager  and  Mack,  1986) . 
Problems  occur  in  defining  the  Laramide  stress  field  in 
southern  New  Mexico  and  Arizona,  as  right-hand  slip  in  New 
Mexico  occurs  on  the  same  fault  systems  as  left-hand  slip  in 
Arizona.  Apparently  there  existed  some  complex  interaction 
between  ENE-WSW  maximum  principal  horizontal  stress  (mphs) 
to  the  southwest  of  the  Colorado  Plateau  (Rehrig  and 
Heidrick,  1976)  and  N-S  mphs  in  the  Rocky  Mountains  (Cries, 
1983)  .  Furthermore,  Drewes  (1981)  and  Bilodeau  (1984) 
report  that  strike-slip  motion  may  involve  pre-existing 
faults,  possibly  faults  associated  with  the  Antler  orogeny 
(Coney,  1978b) .  This  may  also  serve  as  an  explanation  for 
right-lateral  reverse  faults  near  Datil  (Figure  6) .  The 
Hickman  and  Red  Lake  fault  zones  strike  at  a  high  angle  to 


Figure  5.  Map  showing  tectonic  model  for  typical  Laramide 
structures  in  southern  New  Mexico.  The  North 
American  craton  was  driving  southwestward , 
underthrusting  continental  blocks  to  form 
asymmetric  uplifts.  After  Seager  and  Mack  (1986) 
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Figure  6.  Map  showing  late  Laramide  tectonic  configuration, 
including  Eocene  uplifts  and  basins,  eastern 
Arizona  and  western  New  Mexico.  Note  the  high 
angle  between  the  north-northeast  trending  Hickman 
(HFZ)  and  Red  Lake  (RLFZ)  fault  zones  and 
northwest  trending  Love  Ranch  basin.  After  Gather 
and  Johnson  (1986). 
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the  structures  further  south,  but  are  interpreted  as 
Laramide  in  age  (Gather  and  Johnson,  1986) .  But  the  primary 
motion  is  vertical  thrusting,  and  a  significant  amount  of 
crustal  shortening  is  apparent  in  the  region  (Seager  and 
Mack,  1986) . 

The  timing  of  Laramide  volcanism  is  in  some  debate 
within  the  literature.  Coney  (1978a)  suggests  a  complicated 
west  to  east  volcanic  progression,  then  compressional 
deformation,  followed  by  a  retrograde  motion  of  volcanism 
and  plutonism.  Seager  and  Mack  (1986)  suggest  the  presence 
of  a  volcanic  arc  that  persisted  throughout  Laramide  time. 

It  stretched  from  southwestern  New  Mexico  into  east  central 
Arizona.  They  conclude  that  volcanism  postdates  movement  on 
basin-bounding  faults,  but  is  coeval  to  sediments  shed  from 
local  highlands.  Tectonic  activity  resulted  in  a 
progressively  deformed  and  heated  crust  throughout  Laramide 
time.  Coney  (1978a)  infers  that  the  thermal  disturbance  and 
complicated  volcanism  is  the  result  of  the  breakdown  of  a 
underlying  flat  Benioff  zone.  This  breakdown  is  coeval  with 
cessation  of  subduction  of  the  Farallon  plate  at  the  Pacific 
margin. 

The  end  of  Laramide  time  is  generally  regarded  to  be 
synchronous  with  initiation  of  mid-Tertiary  volcanism.  In 
New  Mexico,  Elston  and  others  (1976a)  identify  this  time  at 
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an  angular  unconformity  between  the  Cretaceous  Mesaverde 
group  and  the  Eocene  Baca  formation.  The  earliest  mid- 
Tertiary  volcanic  flows  and  ash-falls  are  interbedded  with 
these  sediments,  but  do  not  affect  Laramide  drainage 
patterns  (Seager  and  Mack,  1986;  Gather  and  Chapin,  1989) . 

The  Baca  formation  has  been  recognized  as  a 
stratigraphic  equivalent  to  the  Eager  formation  and  Mogollon 
Rim  gravels  to  the  west  (Gather  and  Johnson,  1986) .  These 
sediments  were  deposited  in  the  Eocene  Baca  basin,  one  of  a 
series  of  basins  and  uplifts  actively  forming  at  this  time 
in  eastern  Arizona  and  western  and  central  New  Mexico 
(Figure  6) .  The  proximal  location  of  the  Baca  basin  to  the 
Mogollon-Datil  volcanic  field  and  preservation  of  its 
stratigraphic  column  has  proven  beneficial  to  timing  the 
onset  of  mid-Tertiary  volcanism.  The  sediments  were 
deposited  in  a  semiarid  environment  of  alluvial  fans, 
braided  streams,  and  temporary  lakes  in  a  closed  basin 
(Gather  and  Johnson,  1986) .  They  are  deposited 
unconforroably  on  upper  Cretaceous  rocks  eroded  to  a  surface 
of  little  relief,  creating  a  section  up  to  760  m  thick. 
Paleocurrent  directions  indicate  that  the  source  areas  are 
from  surrounding  Laramide  uplifts.  The  Mogollon  Highland  to 
the  south  and  southwest,  the  Sierra  and  Lucero  uplifts  to 
the  east,  and  the  Defiance  and  Zuni  uplifts  to  the  north  all 
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contributed  sediments  to  the  Baca  basin  (Figure  7) .  In  New 
Mexico,  the  sediments  grade  transitionally  into  the  Spears 
formation,  consisting  of  latest  Eocene  andesitic 
volcanoclastic  material.  In  eastern  Arizona,  however,  the 
Mogollon  Rim  gravels  were  eroded,  possibly  beginning  in 
Oligocene  time  (Pierce  and  others,  1979) .  Thus,  the 
stratigraphic  record  in  Arizona  and  New  Mexico  shows  the 
Colorado  Plateau  to  be  topographically  low  with  respect  to 
the  region  to  the  south  until  latest  Eocene. 

Mid-Tertiary 

The  Laramide  orogeny  resulted  in  a  massive 
reorganization  of  the  crust.  This  continued  into  post- 
Laramide  time  with  the  ascent  of  liquids  into  the  upper 
crust,  creating  massive  volcano-plutonic  complexes 
throughout  the  western  United  States  and  Mexico  (Elston  and 
Bornhorst,  1979) .  The  total  volume  of  vented  material  was 
of  sufficient  magnitude  that  structural  controls  for  each 
volcanic  center  were  buried  (Elston  and  others,  1976a) .  In 
all,  more  than  1x10*  km^  of  silicic  material  was  vented  over 
an  area  greater  than  1x10*  km^  in  southwestern  North  America 
(Elston  and  Bornhorst,  1979;  Christiansen,  1989;  Elston  and 
Abitz,  1989). 


Figure  7 .  Map  showing  distribution  of  sediments  off  of  the 
Laramide  Mogollon  highland  and  paleocurrent 
directions  to  the  Mogollon  Rim  gravels,  Eager  fm 
and  Baca  fm.  Southward-directed  paleocurrent 
directions  from  the  Defiance  and  Zuni  uplifts  show 
the  southeastern  Colorado  Plateau  to  be  a 
topographic  low  at  this  time.  After  Gather  and 
Johnson  (1986) . 
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The  surface  features  associated  with  these  complexes 
are  referred  to  in  the  literature  as  ignimbrite  "flareups" 
(Elston,  1976;  Coney,  1978a,  1978b) .  The  exact  cause  of 
these  events  is  unknown,  although  Coney  (1978a)  regards  them 
as  the  result  of  the  retrograde  east-west  sweep  of  igneous 
activity.  These  complexes  developed  during  the  mid- 
Tertiary,  and  in  part  owe  their  existence  to  an  already 
heated  crust  (Coney,  1978a;  Baldridge  and  others,  1989) . 
Silicic  volcanism  occured  as  far  east  as  the  Trans-Pecos 
volcanic  field  during  this  time  (Barker,  1979) .  At  the  same 
time  in  southern  Arizona,  large  metamorphic  core  complexes 
formed  off  the  southern  edge  of  the  Colorado  Plateau. 

Within  this  province,  loci  for  four  mid-Tertiary 
volcano-plutonic  complexes  occur  at  the  margins  of  the 
Colorado  Plateau  (Elston,  1989) .  The  San  Juan,  Colorado 
(Lipman,  1989) ;  White  Mountains,  Arizona;  Marysville,  Utah; 
and  Mogollon-Datil ,  New  Mexico  volcanic  fields  (Elston  and 
others,  1976a)  contain  batholithic  volumes  of  extrusive 
material . 

Of  interest  to  this  study  is  the  Mogollon-Datil 
volcanic  field,  which  lies  within  the  transition  zone 
between  the  Colorado  Plateau,  Basin-Range,  and  Rio  Grande 
Rift  tectonic  provinces.  Previous  geologic  and  geophysical 
work  yields  inconclusive  data  to  make  a  determination  of  the 


relationship  between  this  complex  and  the  plateau.  However, 
this  complex  created  a  volcanic  display  in  the  mid-Tertiary 
that  rivals  any  volcanic  activity  over  the  history  of  the 
planet  (Coney,  197aa) ,  and  its  proximity  to  the  Colorado 
Plateau  is  of  importance  to  unravelling  Cenozoic  history  in 
the  region.  A  more  detailed  review  of  the  Mogollon-Datil 
event  is  provided  later. 

Late  Cenozoic 

The  mid-Tertiary  period  of  ignimbrite  flareup  is 
correlated  in  space  and  time  to  the  transition  from  Laramide 
compression  to  late  Cenozoic  extension.  There  is  little 
debate  that  the  stress  field  in  the  southwestern  United 
States  was  strongly  compressional  during  the  Laramide, 
resulting  in  crustal  shortening  (Seager  and  Mack,  1986) . 
Tension  and  crustal  extension  has  been  dominant  in  the  late 
Cenozoic  throughout  the  region  (Seager  and  others,  1984; 
Henry  and  Price,  1986;  Morgan  and  others,  1986).  Elston  and 
others  (1976a)  state  that  the  high-silica  volcanics  were 
emplaced  in  a  neutral  stress  field. 

The  late  Tertiary  brought  block  faulting,  bimodal 
(mostly  basaltic)  volcanism,  and  extension  in  the 
southeastern  Basin-Range  and  Rio  Grande  Rift  provinces 
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(Coney,  1978a,  1978b;  Morgan  and  others,  1986).  Also, 
crustal  thinning  has  occurred  in  both  provinces  (Warren, 
1969;  Cook  and  others,  1979;  Olsen  and  others,  1979;  Sinno 
and  others,  1981,  1986;  Daggett  and  others,  1986;  Perry  and 
others,  1988) .  Volcanic  fields  are  associated  with  the 
southern  margin  of  the  Colorado  Plateau  in  Arizona  (Luedke 
and  Smith,  1978)  and  the  Jemez  lineament  in  New  Mexico 
(Aldrich  and  Laughlin,  1984) .  The  tectonic  style  associated 
with  this  latest  stage  of  deformation  appears  to  be  distinct 
from  the  previous  ignimbrite  episode,  and  the  chemistry  of 
the  volcanics  indicates  a  different  source  as  well. 

Many  authors  (e.g.  Cather  and  Chapin,  1989)  date  the 
beginning  of  tension  by  the  appearance  of  the  bimodal 
volcanic  suite.  Thus,  tension  preceded  faulting  in  the 
region  (Price  and  Henry,  1984;  Henry  and  Price,  1986). 
East-northeast  extension  apparently  progressed  from  west  to 
east  across  the  region.  Cather  and  Chapin  (1989)  date 
initiation  of  tension  at  36  Ma  in  the  Mogollon-Datil 
volcanic  fie_d,  but  not  until  31  Ma  in  Trans-pecos  Texas 
(Price  and  Henry,  1984;  Henry  and  Price,  1986).  Eocene 
extension  is  supported  by  Cather  and  Johnson  (1986) ,  who 
cite  thickening  of  the  Baca  formation  between  the  Hickman 
and  Red  Lake  fault  zones.  These  faults  are  therefore 
reactivated  Laramide  features,  and  now  bound  a  shallow 
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Eocene  graben.  The  timing  of  reactivation,  however,  is  not 
well-constrained . 

Seager  and  others  (1984)  have  determined  that  two 
periods  of  extension  have  occurred  in  the  Basin-Range  and 
Rio  Grande  Rift  immediately  to  the  southeast  of  the 
Mogollon-Datil  volcanic  field.  The  first  phase  produced 
northwest-trending  moderately  deep  basins  and  uplift  of 
fault-block  mountains.  This  may  have  reflected  a  back-arc 
environment  about  29-28  Ma,  during  a  westward  progression  of 
basaltic-andesite  volcanism. 

The  latter  phase  of  extension,  from  9-3  Ma,  may 
indicate  accelerating  or  renewing  of  block  faulting,  now  in 
a  north-south  trend.  This  is  also  associated  in  time  with 
1-2  km  of  regional  uplift  (Morgan  and  Swanberg,  1985) ,  and 
is  part  of  a  west  to  east  progression  of  block  faulting 
which  began  at  about  20  Ma  in  California.  The  driving  force 
behind  this  phase  of  extension  may  be  strike-slip  tension 
generated  by  the  San  Andreas  transform  zone  to  the  west 
(Seager  and  others,  1984),  and  other  manifestations  include 
lithospheric  thinning,  increased  heat  flow,  and  true 
basaltic  volcanism,  although  true  basalt  is  rare  relative  to 
basaltic  andesite  in  New  Mexico  and  Arizona  (Baldridge  and 


others,  1989) . 
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An  important  point  to  note  is  that  Basin  and  Range 
style  faulting  is  not  present  in  either  the  Colorado  Plateau 
or  the  Mogollon  plateau  (Elston  and  others,  1976a;  Coney, 
1978a) .  This  may  indicate  that  the  Mogollon  plateau  rests 
on  the  Colorado  Plateau,  but  normal  faulting  occurs  between 
the  two  as  well  (Dane  and  Bachman,  1965;  Wilson  and  others, 
1969)  . 

On  the  Colorado  Plateau,  late  Cenozoic  tectonic 
activity  has  been  confined  to  volcanism  near  its  margins 
(Figure  8) .  In  the  study  area,  several  volcanic  centers 
have  developed  along  a  northeast  trending  line  between  the 
White  Mountains-Springerville  and  Jemez  Mountains  volcanic 
fields  (Luedke  and  Smith,  1978) .  This  feature  is  referred 
to  as  the  Jemez  lineament  (Aldrich  and  Laughlin,  1984) . 

With  the  exception  of  the  Zuni-Bandera  field  near  Grants, 

New  Mexico,  Ander  and  Huestis  (1982)  found  no  gravity 
anomalies  associated  with  these  centers.  Furthermore,  Ander 
and  Huestis  presented  magnetotelluric  evidence  that  the 
gravity  "root”  of  the  Zuni-Bandera  field  is  cold,  and  was 
probably  emplaced  before  the  Laramide  orogeny.  Thus,  the 
recent  volcanics  along  the  lineament  are  probably  derived 
from  magma  chambers  beneath  the  upper  crust,  and  are 
composed  of  a  bimodal  suite  of  alkalic  and  tholeiitic 
basalts  (Ander  and  Huestis,  1982) . 


Figure  8.  Map  showing  the  Colorado  Plateau  and  its 

relationship  to  surrounding  provinces.  The 
transition  with  the  extensional  provinces  to  the 
west,  south,  and  southeast  are  marked  by  loci  of 
several  late  Cenozoic  volcanic  fields.  Labelled 
volcanic  fields  include:  S.C.,  San  Carlos;  W.-S. 
White  Mountains-Springerville;  C.C.,  Catron 
County;  Z.-B.,  Zuni-Bandera ;  M.T.-M.C.,  Mount 
Taylor-Mesa  Chevato;  J.M.,  Jemez  Mountains;  T.P., 
Taos  Plateau;  M.C.,  More  County;  and  R.-C.,  Raton 
Clayton.  After  Aldrich  and  Laughlin  (1984) . 


30 


Present-Day  Structure  and  Tectonics 

Colorado  Plateau 

The  Colorado  Plateau  has  remained  relatively 
unaffected  by  tectonic  activity  with  respect  to  the  Basin 
and  Range  province  and  Rio  Grande  Rift  throughout 
Phanerozoic  geologic  time.  The  base  of  the  crust  (Moho) 
under  the  plateau  is  estimated  by  seismic  refraction 
experiments  and  surface  wave  studies  to  be  40-45  km  below 
sea  level  (Roller,  1965;  Toppozada  and  Sanford,  1976;  Keller 
and  others,  1979a) .  However,  crustal  thickness  is  estimated 
to  be  50  km  by  deep  reflection  profiling  (Hauser  and  Lundy, 
1989)  and  earthquake  refraction  data  (Beghoul  and  Barazangi, 
1989)  . 

Upper  mantle  velocities  determined  by  the  seismic 
refraction  surveys  and  surface  wave  dispersion  studies  are 
slow  (7.8  km/sec  at  40  km).  Those  found  at  50  km  (8.1 
km/sec)  are  more  typical  of  the  stable  North  American 
craton.  These  data  indicate  that  the  Moho  under  the  plateau 
may  be  transitional  over  a  10  km  depth  interval,  and  does 
not  appear  as  a  single  reflector.  The  implications  of  this 
possible  transition  are  not  addressed  or  pursued  in  this 
study. 
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The  seismic  refraction  study  of  Jaksha  (1982)  indicates 
a  velocity  of  8.0  km/sec  at  33  km  below  the  Mogollon-Datil 
volcanic  field.  This  result  must  be  viewed  with  caution 
because:  1)  it  has  not  been  determined  whether  the 

Mogo] lon-Datil  field  is  part  of  the  Colorado  Plateau;  and  2) 
with  an  average  of  18  km  station  spacing,  these  data  lack 
the  resolution  of  other  work  in  the  area.  Both  of  these 
points  are  dealt  with  in  the  present  study,  as  this  line  is 
reevaluated  in  light  of  new  seismic  information. 

The  slow  seismic  velocity  at  40  km  indicates  that  the 
upper  mantle  is  presently  hot.  High  electrical 
conductivities  (Pedersen  and  Hermance,  1981)  support  this 
conclusion.  Heat  flow  data  (Reiter  and  others,  1975,  1979; 
Blackwell,  1978;  Thompson  and  Zoback,  1979;  Shearer  and 
Reiter,  1981;  Reiter  and  Clarkson,  1983;  Swanberg  and 
Morgan,  1985)  show  the  interior  of  the  Colorado  Plateau  to 
have  moderate  heat  flow,  about  1.5-1. 8  HFU.  This  is  bounded 
by  a  100  km  wide  zone  of  increasing  heat  flow  towards  the 
Basin  and  Range  and  Rio  Grande  Rift  (Lachenbruch  and  Sass, 
1978;  Swanberg  and  Morgan,  1985).  The  moderate  surface  heat 
flow  does  not  predict  deep  temperature  conditions  implied  by 
seismic  refraction  and  electrical  conductivity  (Morgan  and 
Swanberg,  1985) .  It  is  possible  that  the  base  of  the  crust 
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is  presently  the  site  of  a  thermal  anomaly,  but  the  thermal 
pulse  has  not  yet  reached  the  surface. 

The  present-day  stress  field  of  the  Basin-Range  and  Rio 
Grande  Rift  provinces  also  appears  to  encroach  on  the 
margins  of  the  Colorado  Plateau  (Humphrey  and  Wong,  1983; 
Aldrich  and  Laughlin,  1984) .  Geophysical  studies  at  the 
Great  basin-Colorado  Plateau  boundary  in  Utah  (Keller  and 
others,  1975)  indicate  lithospheric  thinning  and  Basin  and 
Range  signatures  inside  the  plateau  margin.  These  studies 
indicate  that  the  Basin  and  Range  province  is  presently 
growing  at  the  expense  of  the  Colorado  Plateau. 

Despite  this  evidence  for  present  day  tectonism  towards 
the  margin  of  the  Colorado  Plateau,  its  free-air  anomaly  has 
a  near-zero  mean,  indicating  isostatic  equilibrium  (Keller 
and  others,  1979a;  Morgan  and  Swanberg,  1985) .  Inspection 
of  the  Bouguer  gravity  anomaly  and  topographic  elevation 
(Thompson  and  Zoback,  1979)  confirms  that  the  Colorado 
Plateau  must  be  isostatically  compensated,  and  that  the 
compensating  material  must  be  in  the  mantle. 

Southern  Arizona  Basin  and  Range  Province 

Because  of  its  elevation,  southeast  Arizona  has  been 
placed  within  the  Mexican  Highlands  physiographic  sub- 


province  of  the  Basin  and  Range  tectonic  province.  This 
portion  of  the  Basin  and  Range  was  uplifted  in  mid-Tertiary 
time  along  with  the  rest  of  the  region,  but  not  as  much  as 
the  Colorado  Plateau  (Pierce  and  others,  1979) . 

There  is  evidence  that  the  southern  Arizona  portion  of 
the  Basin  and  Range  province  has  unique  characteristics,  and 
therefore  may  be  differentiable  from  the  rest  of  the 
province.  For  example,  the  free  air  gravity  anomaly  map 
shows  a  strong  correlation  with  topography  and  surface 
geology  (Aiken,  1978)  .  This  implies  a  homogeneous  rigid 
lithosphere  in  isostatic  equilibrium.  A  homogeneous  crust 
is  also  inferred  from  heat  flow  data.  Variability  of  all 
heat  flow  readings  in  southern  Arizona  (+  0.46  HFU)  is  less 
than  that  for  the  entire  Basin  and  Range  (+  0.75  HFU).  Heat 
flow  from  wells  greater  than  650  m  in  depth  show  a  southern 
Arizona  heat  flow  of  1.94  +  0.12  HFU  (Shearer  and  Reiter, 
1981),  which  may  be  as  little  as  0.3  HFU  higher  than  the 
Colorado  Plateau  .  Shearer  and  Reiter  state  that  the 
difference  may  be  due  to  the  greater  volume  of  volcanics  and 
plutons  in  the  Basin  and  Range  province,  which  implies  a 
different  lithospheric  response  to  heating  from  below. 

Some  evidence  exists  for  modification  of  the  Basin  and 
Range  crust  in  recent -time.  Sinno  and  others  (1981)  suggest 
that  crustal  thinning  in  southern  Arizona  has  occurred  over 


the  past  5  Ma,  roughly  equivalent  to  the  time  that  north- 
south  structures  were  fully  developed  in  the  southern  Rio 
Grande  Rift  (Seager  and  others,  1984)  and  commencement  of 
regional  uplift  (Morgan  and  Swanberg,  1985) .  Late 
Quaternary  faults  increase  in  frequency  towards  the  east 
side  of  the  Basin  and  Range  into  the  Rio  Grande  Rift, 
implying  that  Basin  and  Range-style  faulting  is  still 
presently  active  in  west  Texas  and  southern  New  Mexico. 

Seismic  velocities  display  evidence  of  elevated  upper 
mantle  temperatures  at  the  present  time.  This  is  based  on 
studies  showing  anomalously  low  Pn  values  throughout 
southern  Arizona;  7.67  km/sec  in  the  southwest  (Sinno  and 
others,  1981)  to  7.82  km/sec  in  the  southeast  (Gish  and 
others,  1981) .  Beneath  the  Mogollon-Datil  volcanic  field, 
near  the  New  Mexico-Arizona  border,  Gish  and  others  (1981) 
found  low  (7.58  km/sec)  Pn  velocities.  As  will  be  discussed 
below,  Pn  velocity  in  the  Rio  Grande  Rift  are  also  low  with 
respect  to  southeastern  Arizona.  The  Basin  and  Range 
province  of  southeastern  Arizona  therefore  has  higher  Pn 
velocities  than  in  the  surrounding  extensional  regions,  but 
still  reflects  values  that  indicate  a  hot  upper  mantle. 

The  residual  Bouguer  gravity  anomaly  of  Arizona  is 
complicated  by  the  overlap  of  shallow  crustal  and  upper 
mantle  heterogeneities  (Aiken,  1978) .  This  also  appears  to 
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be  true  for  the  Basin  and  Range  as  it  enters  southwestern 
New  Mexico  (DeAngelo  and  Keller,  1988) .  In  order  to 
interpret  features  of  a  given  structural  style  or  depth, 
bandpass  filtering  can  be  utilized  (Aiken,  1978) .  However, 
Aiken  also  determined  an  almost  exact  correlation  between  a 
long-wavelength  (800  km)  regional  anomaly  and  the  elevation 
of  gravity  stations  in  southeastern  Arizona.  The  long 
wavelength  regional  trend  was  therefore  removed  by 
cancelling  the  effects  of  topography,  a  procedure  roughly 
equivalent  to  a  high-pass  filter. 

Aiken  (1978)  and  Sumner  (1985)  note  that  residual 
Bouguer  gravity  and  magnetic  trends  are  different  between 
the  Basin  and  Range  and  Colorado  Plateau  provinces.  The 
trends  are  to  the  northwest  in  the  Basin  and  Range,  and  to 
the  northeast  in  the  Colorado  Plateau.  They  are  separated 
by  a  northwest-trending  50  km-wide  transition  zone  in  which 
residual  Bouguer  gravity  anomalies  change  by  approximately 
20  mgal  (Aiken,  1978) .  Residual  anomaly  values  are  in 
general  negative  in  the  Basin  and  Range,  and  positive  in  the 
Colorado  Plateau.  These  patterns  are  thought  to  be 
controlled  by  tectonic  trends  in  the  basement  structure 
(Sumner,  1985) . 

The  physiographic  transition  between  the  Basin  and 
Range  and  Colorado  Plateau  provinces  has  been  placed  at  the 
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Mogollon  Rim  (see  e.g.  Mayer,  1979;  Pierce  and  others, 

1979) .  Based  on  reversed  seismic  refraction  profiling, 
Warren  (1969)  concludes  that  the  Mogollon  Rim  lies  above  an 
abrupt  change  in  crustal  thickness  of  approximately  4  km  in 
central  Arizona.  Additional  thickening  occurs  from 
southwest  to  northeast  as  a  result  of  a  0.5°  dip  on  the 
Moho.  However,  other  tectonic  and  geophysical  indicators 
show  that  the  boundary  is  more  diffuse,  and  that  a 
transition  zone  exists  through  which  geophysical  parameters 
change  (Aiken,  1978;  Brumbaugh,  1987). 

The  geological  and  geophysical  evidence  concerning 
faulting  and  stress,  volcanism,  heat  flow,  gravity  and 
magnetic  trends,  and  seismic  data  all  point  to  the  Basin  and 
Range  province  as  being  a  region  of  tectonic  activity  at 
present.  Many  questions  remain  concerning  this  province. 
Seismic  investigations  are  currently  being  undertaken  along 
the  PACE  (Pacific-Arizona  Crustal  Transect)  in  south-central 
Arizona  (J.  McCarthy,  personal  communication,  1989)  and  from 
Tyrone  through  Lordsburg  to  the  Mexican  border  (S.  Harder, 
personal  communication,  1989) . 
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Rio  Grande  Rift 

The  Rio  Grande  Rift  is  a  complex  region  with  respect  to 
structural  style  and  evolution  (Cook  and  others,  1979; 

Seager  and  Morgan,  1979;  Morgan  and  Golombek,  1984;  Olsen 
and  others,  1987) .  It  is  a  major  continental  rift  system 
(Keller,  1986;  Olsen  and  others,  1987).  However,  it  was 
not  recognized  as  such  until  the  1970s  (e.g.  Olsen  and 
others,  1987;  Keller  and  others,  1989).  At  that  time,  the 
extent  of  the  rift  was  disputed,  with  some  (e.g.  Kelley, 
1979)  extending  it  southward  only  to  about  33°  N  latitude, 
and  others  (e.g.  Gries,  1979;  Seager  and  Morgan,  1979;  Smith 
and  Jones,  1979)  extending  it  as  far  southward  as  northern 
Mexico.  Although  it  corresponds  in  style  and  timing  with 
the  Basin  and  Range  province  in  southern  Arizona,  it  has 
been  shown  to  be  a  separate  feature  (Daggett  and  others, 
1986;  Sinno  and  others,  1986). 

The  region  that  defines  the  Rio  Grande  Rift  is  a  zone 
approximately  1000  km  long  between  central  Colorado  and 
northern  Mexico  (Kelley,  1979) .  It  consists  of  a  series  of 
asymmetric  grabens,  and  forms  the  eastern  margin  of  the 
Colorado  Plateau  in  New  Mexico.  South  of  the  Colorado 
Plateau,  its  surface  expression  widens  to  include  Basin  and 
Range  style  block  faulting  in  southern  New  Mexico,  northern 
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Mexico  and  west  Texas  (Seager  and  Morgan,  1979;  Olsen  and 
others,  1987) . 

The  Rio  Grande  Rift  is  a  region  of  intense  post- 
Laramide  tectonic  activity,  but  is  enigmatic  with  respect  to 
current  tectonic  deformation.  For  example,  late  Quaternary 
faults  frequently  appear  throughout  the  rift  (Seager  and 
Morgan,  1979) ,  yet  contemporary  seismicity  is  not 
differentiable  from  surrounding  provinces  (Jaksha  and 
Sanford,  1986) .  Young  volcanic  centers  with  varying  fluid 
compositions  appear  along  the  rift  throughout  New  Mexico 
(Baldridge,  1979;  Lipman  and  Mehnert,  1979;  Warren  and 
others,  1979;  Zimmerman  and  Kudo,  1979),  and  mid-  and  upper- 
crust  magma  bodies  have  been  identified  (e.g.  Sanford  and 
others,  1977;  Rinehart  and  others,  1979;  Larsen  and  others, 
1986) .  Yet  the  total  volume  of  vented  material  is  minor 
compared  to  the  East  African  Rift  (Olsen  and  others,  1987) . 
Despite  this  seemingly  controversial  evidence,  much 
geophysical  data  exists  that  shows  the  Rio  Grande  Rift  to  be 
an  active  tectonic  feature,  distinct  from  the  provinces  that 
surround  it. 

Quaternary  faulting  and  volcanism  are  related  to  a 
change  in  lithospheric  structure  under  the  Rio  Grande  Rift. 
Many  of  these  changes  are  related  to  heating  at  depth 
(Seager  and  Morgan,  1979).  Heat  flow  readings  exceeding  2.0 
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HFU  are  common.  These  values  are  thought  to  be  caused  by 
transport  of  fluids  along  through-going  fracture  systems  in 
the  crust  (Reiter  and  others,  1979;  Swanberg,  1979;  Reiter 
and  others,  1986) .  Heat  flow  values  are  also  high  near 
volcanic  centers,  yet  are  more  intermediate  away  from  them 
(Reiter  and  others,  1986) .  Due  to  the  lag  time  necessary 
for  conduction  of  heat,  Seager  and  Morgan  (1979)  demonstrate 
that  the  intermediate  background  values  do  not  preclude 
current  heating  at  depth. 

Further  evidence  for  heating  at  depth  comes  from 
seismic  velocity  analysis.  Sinno  and  others  (1986)  show 
inteirmediate  crustal  P-wave  velocities  within  the  rift  when 
compared  to  Basin  and  Range  (BRP)  and  Great  Plains  (GPP) 
provinces.  They  show  upper-middle  crust  velocities  of  5.9- 
6.1  km/sec  in  the  Rio  Grande  Rift  (RGR) ,  as  opposed  to  5.8- 
6.0  km/sec  in  the  BRP  and  6.2  km/sec  in  the  GPP.  Middle 
crust  values  are  6.5  km/sec  in  the  BRP,  6.6  km/sec  in  the 
RGR,  and  6.7  km/sec  in  the  GPP.  They  note  that  lower  crust 
velocity  differences  between  the  BRP  and  RGR  are  near  the 
error  limit  for  their  study,  but  the  change  in  velocity 
could  be  caused  by  injection  of  basalt  dikes  under  the  rift. 
But  upper  mantle  velocities  are  anomalously  low  under  the 
RGR  (7.7  km/sec  as  opposed  to  8.0  km/sec  under  the  adjacent 
BRP  and  8.2  km/sec  under  the  GPP).  This  implies  a 
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significant  amount  of  heat,  and  possible  erosion  of  the 
mantle  lithosphere. 

The  low  Pn  velocity  found  in  seismic  refraction  studies 
for  the  upper  mantle  is  also  reflected  in  S-wave  velocities 
determined  from  surface  wave  dispersion  analyses.  Keller 
and  others  (1979b)  showed  the  Sn  velocity  to  be  about  4.4 
km/sec  (as  opposed  to  4.5  km/sec  in  the  Colorado  Plateau 
(CP)  and  4.6  km/sec  in  the  GPP),  based  on  recordings  from  an 
earthquake  in  Chihuahua,  Mexico.  The  travel  path  for  that 
event,  however,  included  a  portion  of  the  CP.  A  more 
complete  analysis  using  surface  waves  recorded  at  World  Wide 
Standardized  Seismograph  Network  (WWSSN)  stations  at 
Albuquerque  and  El  Paso  (Sinno  and  Keller,  1986)  found  the 
Sn  velocity  to  be  4.25  km/sec  under  the  rift. 

A  local  surface  wave  study  in  the  Albuquerque-Belen 
basins  (Schlue  and  others,  1986)  shows  two  low  velocity 
layers  in  the  crust,  with  the  Moho  interpreted  to  be  close 
to  38  km  depth.  The  19  km  low  velocity  layer  corresponds  to 
the  mid-crustal  magma  chamber  identified  by  Sanford  and 
others  (1977) ,  However,  the  6  km  low  velocity  layer  and 
thickness  of  the  crust  were  unusual  results  for  this  region. 
The  near-surface  low  velocity  layer  was  interpreted  to  be  a 
possible  newly-discovered  magma  chamber.  However,  Schlue 
and  others  (1986)  point  out  that  their  resolvable  data  is 
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restricted  to  periods  of  6-12  sec.  The  analysis  of  Sinno 
and  Keller  (1986)  uses  data  resolvable  between  periods  of  8- 
65  sec,  and  is  in  greater  agreement  with  other  geophysical 
data  for  the  rift  as  a  whole.  But  the  discrepancy  between 
the  two  studies  should  not  be  overly  stressed,  since 
dispersion  between  two  recording  points  reflects  the  average 
velocity  structure  of  the  intervening  material. 

The  studies  of  Keller  and  others  (1979b)  and  Sinno  and 
others  (1986)  also  shows  that  the  crust  is  thinner  under  the 
Rio  Grande  Rift  than  either  side.  The  Moho  in  southern  New 
Mexico  lies  approximately  27  km  below  sea  level,  and  is  at 
about  30  km  and  50  km  under  the  Basin  and  Range  and  Great 
Plains  provinces,  respectively.  Thinned  crust  is  supported 
by  gravity  analysis  (Daggett  and  others,  1986) ,  and  appears 
to  be  thinnest  under  northwest  Chihuahua,  Mexico.  Their 
gravity  modeling  also  indicates  that  the  upper  mantle  under 
the  rift  is  less  dense  than  would  be  expected,  and  they 
interpret  the  cause  to  be  thermal  expansion. 

The  combination  of  thinned  crust,  low-density  upper 
mantle  material,  hot  springs  along  deep  (crustal 
penetrating?)  faults,  low  rigidity  and  slow  upper  mantle 
velocities  all  point  to  the  conclusion  that  heating  is 
occurring  at  the  base  of  the  crust  at  the  present  time.  At 
least  one  stage  of  lithospheric  thinning  occurred  up  to  8  Ma 
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ago  (Perry  and  others,  1988) ,  indicating  that  asthenospheric 
material  has  replaced  mantle  lithosphere  in  the  region 
(Olsen  and  others,  1987;  Perry  and  others,  1987).  The  total 
amount  of  thinning  has  not  yet  been  established.  Based  on 
analysis  of  teleseismic  data,  Davis  and  others  (1989) 
conclude  that  asthenospheric  upwelling  is  presently 
occurring  under  the  Rio  Grande  Rift.  The  amount  of 
upwelling  is  greatest  under  the  southern  part  of  the  rift. 

in  addition  to  crustal  and  lithospheric  thinning,  the 
Rio  Grande  Rift  has  been  elevated  along  with  the  surrounding 
region  (Kelley  and  Duncan,  1986) .  Superimposed  on  this 
uplift  are  areas  of  local  uplift  and  subsidence  (Reilinger 
and  others,  1979) .  For  example,  the  region  near  Socorro  is 
undergoing  active  uplift,  while  downwarp  is  occurring  in  the 
immediate  surrounding  area  (Larsen  and  others,  1986) .  Larsen 
and  others  suggest  that  the  cause  of  this  movement  is  the 
lateral  flow  of  magma  into  the  Socorro  magma  chamber 
combined  with  vertical  fluid  flow  from  depth. 

These  modifications  to  the  lithosphere  beneath  the  Rio 
Grande  Rift  and  resulting  geophysical  signatures  distinguish 
it  from  surrounding  provinces  (Sinno  and  others,  1986;  Olsen 
and  others  1987) .  However,  it  cannot  be  separated  from 
adjacent  regions  within  the  context  of  post-Laramide 
tectonic  events.  For  the  purposes  of  the  present 


discussion,  it  is  important  to  note  the  differences  in 
velocity,  gravity,  crustal  thickness  and  heat  flow  between 
the  Colorado  Plateau  and  the  Rio  Grande  Rift  in  New  Mexico 
It  is  the  transition  between  these  provinces  that  is  the 
target  of  this  study. 


MOGOLLON-DATIL  VOLCANIC  FIELD 


Introduction 

The  Mogollon-Datil  volcanic  field  (MDVF)  is  located  in 
west-southwest  New  Mexico,  near  the  center  of  the  present 
study  (Figure  9) .  It  is  located  in  the  transition  zone 
between  the  Colorado  Plateau,  Basin  and  Range,  and  Rio 
Grande  Rift  provinces.  Some  associated  volcanic  centers 
appear  to  the  south  and  southeast  in  the  extensional 
provinces,  but  none  are  identified  to  the  north  and 
northwest  on  the  plateau. 

The  MDVF  is  part  of  a  province  greater  than  1  x  10^  km^ 
in  area  in  which  Laramide  and  post-Laramide  volcanism  is 
widespread  (Figure  10) .  Elston  (1989)  considers  it  to  be  an 
outlier  of  the  larger  Mexican  Sierra  Madre  Occidental  to  the 
south.  They  are  now  separated  by  the  Mexican  Highlands 
portion  of  the  Basin  and  Range  province  (Figure  10) . 

Volcanism  in  the  MDVF  consists  primarily  of 
intermediate  to  silicic  ash-flow  tuffs  and  lava  flows. 

These  were  vented  by  caldera-type  volcanism  as  part  of  the 
igniinbrite  flareup  in  western  North  America.  At  least  28 
cauldrons  have  been  identified  (Figure  11) ,  and  up  to  35  may 
exist  (Elston,  1978) . 
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Figure  9.  Sketch  map  of  mid-Tertiary  Mogollon-Datil  volcanic 
field  in  southwestern  New  Mexico  and  far  eastern 
Arizona.  After  Ratt6  (1989). 
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Figure  10.  Map  showing  Mogollon-Datil  volcanic  field  as  part 
of  a  mid-Tertiary  volcanic  event  that  involved 
much  of  the  western  United  States  and  Mexico. 
After  Elston  (1978) . 
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Figure  11.  Intermediate  to  silicic  ash-flow  tuff  cauldrons 
associated  with  the  Mogollon-Datil  volcanic 
field,  southwestern  New  Mexico.  For  information 
on  the  numbered  cauldrons,  please  refer  to  Elston 
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General  Geology 

The  MDVF  consists  of  three  volcanic  suites  (Elston  and 
others,  1976a) :  1)  calc-alkalic  andesite  to  rhyolite;  2) 

high  silica  alkali  rhyolite;  and  3)  basalt  and  calc-alkalic 
basaltic  andesite.  These  overlap  in  time  and  space,  yet 
each  appears  to  reflect  different  crustal  and  tectonic 
conditions  at  the  time  of  extrusion  (Figure  12) . 

The  extrusives  of  the  MDVF  are  found  in  and  around  a 
structural  feature  known  as  the  Mogollon  plateau  (MP) , 
described  as  a  giant  resurgent  cauldron  complex  (Elston  and 
others,  1976a) .  In  its  gross  structure,  the  MP  can  be 
divided  into  three  parts:  1)  a  broad,  gently  tilted 
interior  platform  showing  relatively  minor  deformation  with 
respect  to  its  perimeter;  2)  an  encircling  rim  of  uplifted 
material;  and  3)  an  outlying  graben  system  adjacent  to  the 
uplifts  (Coney,  1976a;  Elston  and  others,  1976a) .  The 
volcanic  features  on  the  plateau  have  been  proposed  to  be 
the  surface  expression  of  an  upper  crustal  granitic 
batholith  (Elston  and  others,  1968,  1970,  1976a;  Rhodes, 
1976)  . 

The  volcanic  units  within  the  MP  have  been  collectively 
termed  the  Datil  group  (Dane  and  Bachman,  1965) .  More 
recently,  Gather  and  Chapin  (1989)  have  informally  expanded 


Figure  12.  Timing  and  volume  of  post-Laramide  volcanic 
activity  in  southwestern  New  Mexico  and  its 
relationship  to  tectonic  events  in  the  western 
United  States.  After  Elston  (1976) . 
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the  use  of  the  term  "Datil  group”  to  include  all 
sedimentation  associated  with  mid-Tertiary  volcanism  in 
southwestern  New  Mexico.  Their  "lower"  Datil  group  contains 
sediments  included  in  the  underlying  Eocene  Baca  formation. 
The  stratigraphic  record  indicates  a  smooth  transition  from 
the  Baca  formation,  and  the  volcanics  are  grouped  as 
andesite  (with  high  KjO) .  The  "upper"  Datil  group  is 
composed  entirely  of  volcanic  deposits,  now  of  bimodal 
composition.  The  two  rock  types  representative  of  this 
group  are  mafic  lavas,  termed  basaltic  andesites,  and 
silicic  ash-flow  tuffs.  The  transition  is  dated  at  36  Ma 
and  is  thought  to  correspond  to  the  change  between  Laramide 
extension  and  Cenozoic  tension  in  the  region. 

Timing  of  Volcanism 

Volcanism  was  initiated  in  mid-Tertiary  time  with  the 
extrusion  of  the  calc-alkalic  andesite  to  rhyolite  group 
(Elston  and  others,  1976a) .  Minor  amounts  of  this  suite 
began  to  vent  at  43  Ma,  with  a  significant  increase  in 
volume  at  36  Ma  (McIntosh,  1989) .  Source  fluids  for  these 
rocks  may  have  been  derived  from  the  remnants  of  a  subducted 
slab,  with  probable  upper  mantle  and/or  lower  crustal 
contamination.  Elston  and  others  (1976a)  consider  this 
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suite  to  be  representative  of  Colorado  Plateau  margin-type 
volcanic  material. 

The  last  group  of  rocks  to  appear  were  the  basalt  and 
calc-alkaline  basaltic  andesites,  erupted  from  37-0  Ma. 

These  are  associated  by  Elston  and  others  (1976a)  to  be 
typical  of  Basin-Range  volcanism.  They  are  thought  to  be  - ^ 
mantle  origin,  caused  by  the  relaxation  of  the  western  North 
American  plate  after  subduction  ceased  on  its  western  margin 
(Elston,  1976;  Elston  and  others,  1976a) .  Asthenospheric 
material  was  allowed  to  move  into  the  lithosphere  at 
liquidus  temperatures  during  this  time.  Thus  these  groups 
are  of  different  origin,  although  they  intersect  in  space 
and  time. 

Overlapping  both  of  these  suites  in  time  (32-21  Ma)  are 
the  rocks  of  the  high  silica  alkalic  rhyolite  suite  (Elston, 
1976;  Elston  and  others,  1976a;  McIntosh,  1989) .  These  form 
a  distinctive  group  of  rocks  with  respect  to  chemistry, 
location  of  volcanism,  and  possibly  source  location  and 
composition.  With  minor  exceptions,  they  occur  entirely 
within  the  Mogollon  plateau. 

The  bulk  of  magmatic  activity  occurred  in  late  Eocene- 
early  Oligocene  time.  McIntosh  (1989)  uses  ^°Ar/^’Ar 
chronology  to  show  four  major  volcanic  episodes,  each 
separated  by  a  hiatus,  from  36.1  to  26.3  Ma.  In  general. 
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volcanic  activity  progressed  from  east  to  west  and  south  to 
north.  Three  of  the  four  episodes  correlate  in  time  to 
volcanism  in  the  San  Juan  volcanic  field  in  Colorado, 
implying  consistent  tectonic  control  on  the  east  side  of  the 
mid-Tertiary  Colorado  Plateau. 

Mogollon  Plateau 

With  the  exception  of  about  2.5  km^  of  exposed 
underlying  material,  the  MP  is  covered  by  volcanic  rocks  of 
Oligocene-Pliocene  age  (Coney,  1976a) .  The  eruptive  style 
for  this  structure  is  that  of  a  series  of  overlapping  flows, 
ash-falls,  and  pyroclastic  and  epiclastic  debris  (Elston  and 
others,  1976a;  Rhodes,  1976).  The  total  thickness  of  the 
volcanic  pile  averages  3.1  km  (Coney,  1976a;  Rhodes,  1976). 
These  units  originate  from  overlapping  resurgent  cauldrons, 
ring-dike  complexes  and  rhyolite  domes. 

Volcanic  deposits  within  the  MP  unconformably  overlie 
Precambrian  crystalline  rocks,  Paleozoic-Mesozoic 
sedimentary  rocks,  and  volcanic  rocks  of  late  Cretaceous  to 
early  Tertiary  age  (Krohn,  1976) .  The  bulk  of  the  volcanics 
on  the  MP  is  composed  of  high  silica  alkali  rhyolite  (Elston 
and  others,  1976a) . 
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There  are  several  noteworthy  geological  and  geophysical 
aspects  of  the  Mogollon  plateau: 

1)  With  the  exception  of  two  minor  centers  of  high 
silica  extrusives  (Red  Mountain,  Arizona,  and  the  Mule  Creek 
caldera.  New  Mexico) ,  all  of  the  high  silica  rocks  are 
clustered  in  a  region  to  the  south  of  the  Plains  of  St. 
Augustine  and  to  the  north  of  silver  City  (Elston  and 
others,  1976a) . 

2)  There  appear  to  be  two  types  of  felsic  rocks  within 
this  suite:  (a)  A  widespread,  phenocryst-poor  "framework 
lava";  and  (b)  A  suite  of  phenocryst-rich  lavas  with  varying 
composition  that  is  spatially  related  to  a  source  cauldron- 
"cauldron  lava"  (Rhodes,  1976) .  Sources  for  the  framework 
lavas  become  younger  as  the  center  of  the  MP  is  approached, 
indicating  a  common,  shrinking  source  region  (Elston  and 
others,  1976a) . 

3)  During  the  time  of  high  silica  volcanism,  the  calc- 
alkalic  and  basalt  to  calc-alkalic  basaltic  andesite  suites 
were  not  vented  on  the  MP.  Since  they  predate  and  postdate 
the  rhyolite  volcanism  on  the  plateau,  Elston  and  others 
(1976a)  suggest  that  a  geochemical  barrier  and/or  thermal 
barrier  were  present  under  the  MP  from  32-21  Ma. 

4)  Gravity  studies  (Coney,  1976a;  Krohn,  1976)  show  low 
values  over  the  MP  that  cannot  be  completely  explained  by 
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surface  geology.  They  interpret  the  anomaly  to  represent  a 
mass  deficiency  in  the  upper  crust. 

5)  The  exterior  grabens  surrounding  the  MP  are  of 
Basin-Range  type  (Elston  and  others,  1976a) .  However,  the 
grabens  do  not  strike  north-northwest  as  do  other  Basin- 
Range  faults  in  the  region.  They  appear  to  wrap  around  the 
MP.  Elston  and  others  (1976a)  interpret  this  to  be  an 
indication  of  some  type  of  barrier  to  faulting. 

Taking  these  observations  together,  Elston  and  others 
(1976a)  have  deduced  that  a  shallow  granitic  plutonic 
complex  exists  under  the  MP.  Rhodes  (1976)  states  that  the 
framework  lavas  represent  liquid  from  the  main  body  of  an 
individual  batholith,  but  cauldrons  represent  areas  above 
solitary,  volatile-rich  cupolas. 


DATA  COLLECTION  AND  COMPILATION 


Seismic  Refraction 

Two  previous  seismic  refraction  studies  were  conducted 
across  the  Mogollon-Datil  volcanic  field  (Figure  13) . 

Jaksha  (1982)  reported  on  a  partially  reversed  east-west 
profile  between  Morenci,  Arizona  and  the  White  Sands  Missle 
range,  New  Mexico.  Harden  (1982)  collected  data  from  an 
unreversed,  north-northeast  trending  line  between  Tyrone  and 
the  Acoma  Pueblo,  New  Mexico.  Information  from  both  of 
these  studies  were  compiled  for  reanalysis  in  this  study. 

Between  July  1987  and  August  1989,  the  Harden  (1982) 
line  was  reoccupied.  The  total  number  of  stations  deployed 
along  this  profile  constituted  a  twofold  increase  of 
recorded  data  over  the  previous  study,  and  an  increase  of 
approximately  125%  over  the  MDVF  as  a  whole. 

Sources  for  this  survey  were  quarry  blasts  at  the 
Phelps  Dodge  open-pit  copper  mine  at  Tyrone,  New  Mexico. 
Since  it  was  not  possible  to  determine  exact  origin  times 
from  blast  tickets  provided  by  Phelps  Dodge,  a  base  station 
was  installed  on  the  mine  property  near  the  pit.  Distance 
from  the  blasts  to  this  base  station  never  exceeded  3  km. 
Velocity  of  P-waves  within  the  pit  was  determined  to  be  5.5 
3cm/sec  (S.  H.  Harder,  personal  communication,  1987) 
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Figure  13.  Location  of  seismic  refraction  lines  used  in  the 
present  study.  The  east-west  trending  line  is 
from  Jaksha  (1982)  .  Locations  of  seismic 
recording  stations  are  indicated  by  a  (+) . 
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from  a  local  study.  With  two  exceptions,  the  base  station 
was  equipped  with  Sprengnether  DR-100  recording  instruments, 
with  L4  geophones  of  1  Hz  natural  frequency.  These  records 
were  stored  on  cassette  tape  and  later  replayed  with  a  DP- 
100  playback  unit.  The  exceptions  were  two  instances  in 
which  portable  Sprengnether  MEQ-800  analog  recorders  were 
used,  records  were  stored  on  smoked  paper  and  later  visually 
inspected.  Arrival  time  errors  were  never  greater  than  .1 
sec  (smoked  paper) ,  and  usually  negligible  (tape  records) . 
Source  parameters  are  provided  in  Appendix  I. 

Stations  were  occupied  by  Sprengnether  MEQ-800  portable 
seismograph  stations.  Portable  FM  tape  recorders  were 
connected  to  the  MEQs  to  record  the  signal.  Sprengnether  S- 
7000  and  L4  seismometers,  both  of  1  Hz  natural  frequency, 
were  used  in  this  study.  Detailed  data  describing  the 
stations  are  presented  in  Appendixes  II  and  III. 

In  all,  59  distinct,  good  quality  seismograms  were 
recorded.  Several  sites  were  reoccupied  in  order  to  insure 
accuracy  in  source-receiver  travel  time  and  in  order  to 
insure  good  signal-to-noise  ratios.  In  all  cases  where 
possible,  the  FM  tapes  were  replayed  into  a  DR-100  for 
digital  sampling  of  the  signal  and  then  input  as  data  files 
on  the  UTEP  Department  of  Geological  Sciences  Hewlett- 
Packard  HP-9000  computer.  When  data  was  unretrievable  from 
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tape  via  this  method,  the  signal  was  hand  digitized,  then 
resampled  using  a  cubic  spline  routine.  Signals  used  in 
this  analysis  were  digitized  at  either  .01  or  .02  sec.  An 
explanation  of  procedures  for  data  retrieval  is  provided  in 
Appendix  IV.  Each  record  was  normalized  to  its  greatest 
amplitude  for  display,  so  amplitudes  of  phases  are  not 
comparable  from  trace  to  trace.  However,  relative 
amplitudes  are  comparable  within  each  trace  in  this  study. 

These  records  were  then  merged  with  the  data  from 
Harden  (1982)  to  produce  the  record  section  shown  in  Figure 
14,  and  the  filtered  record  section  in  Figure  15.  Arrivals 
of  significant,  recognizable  phases  were  picked  for  each 
record  and  input  to  a  ray  tracing  program  for  analysis. 

Data  from  the  Jaksha  (1982)  cross-line  (Figure  16)  was 
used  in  the  present  study  for  better  control  on  three- 
dimensional  structure.  Arrivals  were  repicked  from  his 
record  sectio.i,  at  an  estimated  error  of  less  than  0.2  sec. 

A  reversed  profile  for  Pg  and  Pn  from  some  of  the  stations 
was  included  in  that  study.  The  Pn  branch  was  repicked  with 
a  similar  estimated  error  as  above,  and  used  to  constrain 
the  upper  mantle  velocity  and  location  of  the  Moho  in  this 
study . 
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Figure  16.  Seismic  refraction  line  of  Jaksha  (1982), 
original  picks  used  in  that  study. 


showing 
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Gravity 

The  extensive  gravity  data  base  at  the  University  of 
Texas  at  El  Paso  is  well-documented  (Daggett  and  others, 
1986;  Jenkins  and  Keller,  1989).  This  information  source 
has  been  utilized  in  two  ways:  1)  to  make  gravity  maps  of 
the  study  area;  and  2)  to  extract  gravity  profiles  along 
each  of  the  seismic  refraction  lines  described  above. 

Complete  Bouguer  gravity  anomaly  values  were  retrieved 
from  a  region  enclosed  by  31.0°  to  36.5°  N  latitude  and 
105.5°  to  112.5°  W  longitude.  A  sea  level  datum  and  a 
density  of  2.67  gm  cm'^  were  used  in  the  Bouguer  and  terrain 
corrections  for  this  study.  The  total  number  of  available 
readings  exceeds  40,000  in  this  area.  Interpolated  gravity 
values  were  derived  at  evenly-spaced  (4  km)  locations  in  the 
region,  using  the  technique  of  minimum  curvature  (Briggs, 
1974;  Swain,  1976).  To  insure  stability  in  frequency-domain 
filtering,  a  third-order  polynomial  surface  was  fitted  to, 
and  subsequently  subtracted  from  the  data  (Figure  17) .  To 
separate  trends  for  analysis,  these  data  were  then  filtered 
using  the  technique  of  Peeples  and  others  (1986) . 

Broad  trends  reflecting  Moho-related  features  required 
a  125  km  low-pass  filter  in  the  wave  number  domain  in  an 
area  overlapping  that  of  the  present  study  (Daggett  and 
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others,  1986) .  This  filter  has  been  applied  and  the  results 
are  shown  in  Figure  18.  Upper  crustal  features  in  the  study 
area  are  identified  using  a  15-125  km  band  pass  filter 
(Figure  19) . 

Gravity  profiles  were  extracted  from  the  data  base 
along  strike  of  the  seismic  refraction  lines  used  in  this 
study.  These  consist  of  Bouguer  anomaly  values  projected 
from  gravity  stations  no  more  than  1  km  from  the  profile. 
They  are  used  in  conjunction  with  the  refraction  lines  for 
interpretation  of  the  deep  structure  of  the  Mogollon-Datil 
volcanic  field.  Plots  of  the  profiles  are  presented  in  the 
discussion  section  which  follows. 

Heat  Flow 

An  extensive  national  heat  flow  data  base  has  been 
provided  for  the  present  study  (P.  Morgan  and  J.  Sass, 
personal  communication,  1988) .  Data  for  Arizona  and  New 
Mexico  have  been  extracted  and  are  shown  in  Figure  20.  This 
figure  supports  previous  work  which  shows  that:  1)  the  heat 
flow  within  the  interior  of  the  Colorado  Plateau  is 
relatively  normal  with  respect  to  stable  continental  values; 
and  2)  the  heat  flow  in  the  Rio  Grande  Rift  and  Basin  and 
Range  provinces  are  above  normal,  with  highest  values 
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occurring  near  the  boundary  of  the  Colorado  Plateau;  3) 
within  these  provinces,  heat  flow  is  highest  near  deep 
faults  and  late  Cenozoic  volcanism;  and  4)  groundwater 
circulation  has  complicated  heat  flow  analysis  in  the 
region.  These  points  are  discussed  in  greater  detail  later. 
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Figure  17.  Bouguer  gravity  roap  of  eastern  Arizona  and 

western  New  Mexico.  A  third-order  polynomial 
surface  fitted  to  the  gridded  values  has  been 
removed  to  produce  this  surface.  Contour 
interval  =  5  mgal. 
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Figure  18.  Regional  Bouguer  gravity  anomaly  map  of  eastern 
Arizona  and  western  New  Mexico  constructed  by  a 
125  km  low-pass  filter  in  the  wave  number  domain. 
Contour  interval  =  5  mgal. 
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Figure  19 .  Regional  Bouguer  gravity  anomaly  map  of  eastern 
Arizona  and  western  New  Mexico  constructed  by  a 
15-125  km  band-pass  filter.  Contour  interval  =  5 
mgal . 
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Figure  20.  Heat  flow  readings  in  Arizona  and  New  Mexico. 

oVTnbols  indicate  the  corresponding  range  of  heat 
flow  values:  pluses  indicate  <  0.8  HFU; 
triangles  correspond  to  0.9  to  1.6  HFU;  hexagons 
indicate  1.6  to  2.4  HFU;  and  diamonds  correspond 
to  values  >  2.4  HFU. 


DATA  ANALYSES,  MODELING  AND  DISCUSSION 


Seismic  Refraction  and  Gravity  Profiling 


Technique 

The  Phelps-Dodge  copper  mine  at  Tyrone  has  been  the 
source  for  several  seismic  refraction  studies  in 
southwestern  New  Mexico  and  southeastern  Arizona  (Figure 
21) .  The  crustal  depth  and  velocity  structure  reported  by 
Sinno  and  others  (1986)  was  used  as  a  constraint  near  Tyrone 
for  the  present  study.  The  intersection  of  the  Jaksha 
(1982)  refraction  line  with  this  profile  has  proven 
beneficial  in  three  ways:  1)  resulting  models  must  be 
consistent  at  the  point  of  intersection,  providing  an 
additional  constraint  to  layer  thickness  and  velocity  in 
each;  2)  the  ambiguity  due  to  the  lack  of  reversal  of  the  N- 
S  profile  is  reduced;  and  3)  they  cross  within  the  Mogollon 
plateau,  a  major  target  of  this  study. 

Gravity  profiles  have  been  chosen  to  closely  follow  the 
seismic  refraction  lines  described  above.  For  each  profile, 
gravity  readings  within  1  km  were  included  by  projection 
along  a  line  perpendicular  to  its  axis.  Thus  gravity  data 
are  used  as  an  independent  constraint  to  seismic  model 
parameters.  The  objective  of  this  analysis  is  to  create  an 
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Figure  21.  Location  map  for  seismic  analyses  of  the  tectonic 
provinces  and  their  boundaries  in  Arizona  and  New 
Mexico.  Sources  for  refraction  data:  a  (this 
study);  b  (McCullar,  1977;  Cook  and  others, 

1978) ;  c  (Toppozada  and  Sanford,  1976) ;  d  (Olsen 
and  others,  1979);  e  (Jaksha,  1982);  f  (Stewart 
and  Pakiser,  1962);  g  (Roller,  1965);  h  (Gish  and 
others,  1981) ;  i  (Sinno  and  others,  1986) . 
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earth  model  satisfying  both  of  these  data  sets,  the  desired 
result  being  an  image  of  present-day  crustal  structure. 

The  seismic  modeling  program  used  to  determine  earth 
structure  is  based  on  the  algorithm  described  by  Cerveny  and 
Ravindra  (1971)  .  This  technique,  as  applied  by  Leutgert 
(personal  communication,  1989;  see  Hill  and  others,  1985) , 
generates  a  model  consisting  of  two-dimensionally  varying 
earth  structures.  Many  workers  have  used  this  ray-tracing 
technique,  called  asymptotic  ray  theory  (ART) ,  and  it  is 
well-described  in  the  literature. 

The  computer  program  used  to  interpret  the  gravity 
profiles  across  the  area  is  based  on  a  2.5  dimensional 
algorithm.  The  program  is  based  on  the  technique  of  Talwani 
and  others  (1959)  and  Cady  (1980).  It  was  written  by  S.  F. 
Lai  at  the  University  of  Texas  at  Dallas,  and  modified 
extensively  by  D.  Roberts  (personal  communication,  1989) . 

Depth  to  major  velocity/density  interfaces  were  the 
same  for  the  seismic  and  gravity  models  at  the  point  of 
intersection  of  the  Tyrone-Acoma  and  Morenci-Dice  Throw 
profiles,  except  for  the  near  surface-upper  crustal 
boundary.  This  exception  is  for  the  following  reasons: 

1)  The  seismic  ray-tracing  program  is  capable  of 
including  topography,  while  gravity  data 
processing  assumes  reduction  to  a  datum  plane 
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effectively  at  the  lowest  station  elevation.  The 
Bouguer  correction  results  in  a  D.C.  shift  equal 
to  the  portion  of  the  correction  due  to  material 
between  this  elevation  and  sea  level. 

2)  Gravity  and  seismic  stations  do  not  coincide. 
Therefore,  the  gravity  readings  may  reflect  very 
local  structure  that  does  not  influence  the 
seismic  ray  path. 

3)  Near-surface  inhomogeneities  within  the 
Mogollon  plateau  may  cause  minor  disagreements 
between  the  two  data  sets. 

The  near  surface  structure  is  smoother  in  the  seismic 
model  than  in  the  gravity  model  because  of  limitations  in 
the  ART  approach.  Because  of  computer  time  limitations,  it 
was  not  possible  nor  desirable  to  model  highly  complex 
shallow  structures.  The  near  surface  structure  in  the 
gravity  models  is  more  complex  than  required  by  either  the 
gravity  or  seismic  data,  but  this  complexity  reflects  known 
geologic  features.  Thus,  the  gravity  models  should  be 
considered  large  scale  geologic  cross-sections. 

Known  geologic  information  was  used  to  create  the 
initial  model  and  constrain  its  evolution.  When  alterations 
to  structure  were  suggested  for  a  particular  profile 
(gravity  or  seismic,  Tyrone-Acoma  or  Morenci-Dice  Throw) 
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during  modeling,  the  effects  on  the  other  profiles  were 
examined.  Thus  every  modification  required  running  two  ART 
and  two  Talwani  programs  to  observe  its  effect.  More  than 
250  repetitions  of  these  computer  programs  were  run  before  a 
satisfactory  earth  structure  was  derived.  The  final  model 
satisfies  both  refraction  and  gravity  observations,  ties  to 
nearby  results  (Olsen  and  others,  1979;  Gish  and  others, 
1981;  Sinno  and  others,  1986)  and  is  compatible  with  surface 
geologic  features. 

Refraction  Observations 

Observed  arrivals  for  major  crustal  velocity  interfaces 
are  shown  respectively  for  the  data  of  Jaksha  (1982)  and  the 
present  study  in  Figures  22  and  23.  For  this  discussion, 
the  following  notations  are  used  to  identify  different 
phases:  Pg  and  Pn  are  critically  refracted  P  waves  from  the 

top  of  the  upper  crust  and  upper  mantle  respectively;  and 
PgR  and  PcR  indicate  reflections  from  the  top  of  the  upper 
and  lower  crust  respectively.  PmP  corresponds  to  the  P-wave 
reflection  from  the  M-discontinuity . 

The  Jaksha  profile  (Figure  22a)  shows  Pg  to  be  a  first 
arrival  to  a  distance  of  155  km,  after  which  the  Pn  branch 


Figure  22.  a)  Repicked  arrivals  for  Jaksha  (1982)  Morenci- 
Dice  Throw  refraction  line, 
b)  Arrival  picks  for  reversed  Pn  branch. 
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Figure  23.  Arrival  picks  for  Tyrone-Acoma  refraction  line. 
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becomes  the  first  arrival.  PgR,  PcR  and  PmP  never  become 
first  arrivals  on  this  line,  yet  their  arrival  times 
combined  with  gravity  data  and  the  intersection  with  the 
Tyrone-Acoma  line  (this  study)  are  important  contributors  to 
determining  crustal  structure.  Upper  mantle  velocity 
structure  and  placement  of  the  Moho  is  supported  by  the  use 
of  the  reversed  Pn  branch  (Figure  22b) . 

The  density  of  seismic  readings  along  the  Tyrone-Acoma 
profile  gives  confidence  to  location  and  timing  of  arrivals. 
Analysis  of  Figure  23  shows  several  items  of  interest 
throughout  the  profile.  Pg  arrivals  are  delayed  in  the 
interval  from  30  km  to  85  km,  indicating  slow  upper  crustal 
velocities.  Chatter  in  the  Pg  data  is  interpreted  to  be 
caused  by  topography  and  near-surface  inhomogeneities  in  the 
volcanic  pile  within  the  southern  Mogollon  plateau. 

PgR  does  not  appear  as  a  distinct  arrival  on  the 
Tyrone-Acoma  line.  Possible  reasons  for  this  include:  1) 

The  upper  crust-middle  crust  interface  does  not  produce  a 
strong  reflector;  2)  The  delay  in  Pg  arrivals  may  have 
interfered  with  the  identification  of  PgR;  and  3)  The  time 
shifts  due  to  near  surface  effects  make  reliable  phase 
correlation  impossible.  PcR  never  occurs  as  a  first  arrival 
in  this  region,  yet  it  gives  valuable  information  that 
provides  a  constraint  to  the  gravity  model  derived  for  this 
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profile.  Pn  becomes  the  first  arrival  at  approximately  170 
km,  and  ceases  to  be  observed  beyond  220  km. 

The  most  striking  feature  in  the  Tyrone-Acoma  profile 
is  the  time  shift  in  PmP  arrivals  from  210  km  to  225  km. 
Across  this  interval,  PmP  arrivals  step  almost  1.5  sec 
forward  on  the  reduced  (T  -  X/6  sec)  travel  time  plot. 

Early  arrivals  indicate  that:  1)  seismic  rays  now  pass 
through  faster  material;  and/or  2)  the  distance  traveled 
has  been  shortened.  A  change  of  fundamental  importance 
within  the  crust  has  created  this  step,  and  must  be 
accounted  for  in  the  final  integrated  model. 

Gravity  Observations 

Gravity  profiles  for  both  refraction  lines  are  shown  in 
Figure  24.  A  gravity  low  is  apparent  in  the  center  of  the 
Jaksha  line  (Figure  24a) .  It  is  approximately  80  mgal  lower 
than  the  west  end  and  70  mgal  lower  than  the  east  end 
(disregarding  the  Rio  Grande  graben  anomaly  at  270  km 
distance) .  By  itself,  this  low  does  not  constrain  the 
density  difference  or  depth  of  the  causitive  anomalous  body. 
However,  the  half  wavelength  (about  250  km)  suggests  a  deep, 
possibly  Moho-related  feature.  A  depression  or  sag  is 
observed  within  this  low  from  about  50  km  to  165  km.  This 
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feature,  about  115  kiti  in  half  wavelength,  may  be  caused  by 
an  upper  crustal  anomalous  body.  The  Tyrone-Acoma  gravity 
field  also  displays  a  low  in  the  center  of  the  profile 
(Figure  24b) .  Values  in  the  center  of  this  feature  are 
nearly  100  mgal  lower  than  those  at  the  south  end,  and 
nearly  50  mgal  lower  than  the  north  end.  Thus,  it  is 
superimposed  on  a  slope  in  the  regional  gravity  field  that 
becomes  more  negative  to  the  north.  Along  the  Morenci-Dice 
Throw  line,  a  subtle  depression  is  observed  from  50  km  to 
200  km  within  the  larger  negative  anomaly. 

Model  Generation  and  Evolution 

Jaksha  (1982)  used  a  flat-lying,  homogeneous  layered 
model  to  interpret  his  results  for  the  Morenci-Dice  Throw 
profile.  His  model  was  approximately  duplicated  using  ART 
modeling  and  the  result  is  shown  in  Figure  25.  Although 
this  structure  is  adequate  in  reproducing  observed  arrivals, 
it  yields  a  Bouguer  anomaly  field  that  does  not  match 
observed  values  (Figure  26) .  Since  the  gravity  field 
indicates  greater  complexity  within  the  crust  and  upper 
mantle,  alternative  models  were  generated. 

The  first  alternative  is  one  in  which  the  gravity 
anomaly  is  derived  entirely  by  crustal  thickening  (Figure 
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27) .  This  model  shows  that  a  near-surface  pluton  is  not 
necessary  to  satisfy  gravity  observations.  The  resulting 
model  was  then  used  in  the  ART  seismic  ray-tracing  program, 
and  the  result  is  shown  in  Figure  28. 

For  the  Morenci-Dice  Throw  profile,  it  appears  from 
seismic  and  gravity  data  that  crustal  thickening  is  a  viable 
explanation  of  geophysical  anomalies  in  the  region.  The 
maximum  thickness  of  the  crust  does  not  exceed  37  km  along 
this  line,  which  is  intermediate  between  the  Colorado 
Plateau  and  Basin  and  Range/Rio  Grande  Rift  provinces.  The 
effect  of  the  upper  mantle  density  change  under  the  Rio 
Grande  Rift  is  to  slightly  thin  the  crust  at  that  point. 
However,  crustal  thinning  would  still  be  observed  even  if 
the  upper  mantle  density  were  made  homogeneous,  and  the 
outcome  of  the  analysis  is  unchanged. 

The  margins  of  crustal  thickening  correspond 
approximately  to  the  -20  mgal  contour  line  of  the  low-pass 
gravity  map  in  Figure  18.  If  this  association  holds  in 
southwestern  New  Mexico,  the  thickened  crust  is  present 
nearly  as  far  south  as  Silver  City.  Therefore,  a  model 
based  purely  on  thickened  crust  indicates  that  the  region 
between  Morenci  and  Dice  Throw  is  part  of  the  transition 
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Figure  25.  Ray  trace  model  using  Jaksha's  (1982) 

interpretation  of  crustal  structure  from  the 
Morenci-Dice  Throw  refraction  line.  Velocities 
are  in  km/sec.  Symbols:  o  =  observed  arrivals:  + 
=  calculated  reflection  arrivals;  x  =  calculated 
refraction  arrivals. 
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Figure  26.  Gravity  profile  based  on  Jaksha's  (1982)  earth 
model.  Pluses  (+)  indicate  observed  values; 
triangles  reflect  model-calculated  gravity  field. 
The  flat-lying  structure  obviously  does  not  fit 
the  observed  data.  Values  for  density  are  in 
gm/cm^. 
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Figure  27.  Gravity  model  for  thickened  crust  showing  the  fit 
of  predicted  values  to  observed  data  along  the 
Jaksha  (1982)  refraction  line.  Pluses  (+) 
indicate  observed  values;  triangles  reflect 
model-calculated  gravity  field.  Crustal 
thickening  is  sufficient  to  account  for  the 
observed  anomaly.  Values  for  density  are  in 
gm/cm^ . 
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The  second  alternative  investigates  whether  a  near¬ 
surface  pluton  could  account  for  the  gravity  anomaly  instead 
of  crustal  thickening.  Of  primary  concern  for  this  model  is 
the  tradeoff  between  density  contrast  and  volume.  The 
density  of  a  granitic  pluton  is  almost  always  lower  than 
that  of  the  displaced  country  rock  (Bott  and  Smithson, 

1967),  usually  within  the  range  of  -0.05  to  -0.18  gm/cm^. 
Bott  and  Smithson  find  that  the  departure  is  rarely  more 
than  -0.2  gm/cm^.  For  the  present  study,  the  density 
contrast  of  the  pluton  is  set  equal  to  -0.1  gm/cm^  with 
respect  to  country  rock  (in  a  lateral  sense) ,  a  conservative 
value  within  their  range. 

The  volume  of  the  plutonic  complex  must  now  be 
considered.  E.  A.  Anthony  (personal  communication,  1989) 
states  that  venting  of  a  pluton  into  a  caldera  volcanic 
complex  may  remove  from  10-90%  of  the  original  volume  of  the 
magma.  If  the  plutonic  complex  of  Rhodes  (1976)  and  Elston 
and  others  (1976a)  covers  the  same  lateral  extent  as  the 
outline  of  the  caldera  complex,  the  thickness  of  the 
underlying  magma  chamber  can  be  calculated  using  the  formula 
for  the  gravitational  attraction  of  a  cylinder.  Given  an 
average  thickness  of  3.1  km  of  volcanic  pile,  90%  removal  of 
liquid  from  the  underlying  chamber  yields  an  approximately 
0.33  km  thick  pluton.  This  pluton  is  too  thin  to  have 
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significant  effect  on  the  gravity  field  or  the  wide  angle 
seismic  measurements.  Venting  of  10%  of  the  magma  chamber 
yields  a  residual  pluton  about  27.9  km  thick.  Since  this 
represents  almost  the  entire  thickness  of  the  crust,  the  10% 
scenario  was  also  discarded  for  this  study. 

Two  models  were  therefore  developed  with  plutons  of 
intermediate  thicknesses:  1)  50%  of  the  magma  was  vented  to 
the  surface,  leaving  a  3.1  km  thick  pluton  in  the 
subsurface;  and  2)  20%  of  the  magma  was  removed,  yielding  a 
12.5  km  thick  pluton.  The  top  of  each  pluton  was  set  at  4.5 
km  below  the  present  surface,  primarily  to  leave  room  for 
the  volcanic  pile  that  presently  exists  in  the  area. 

Figure  29  displays  the  3.1  km  pluton  in  the  upper 
crust,  and  shows  that  the  gravity  anomaly  produced  is 
insufficient  for  observed  readings.  This  model  would 
require  either  a  greater  density  contrast  with  respect  to 
the  country  rock,  or  a  component  of  crustal  thickening  to 
supply  mass  deficiency  with  depth.  Figure  30  shows  that  the 
mass  deficiency  of  a  pluton  about  12.5  km  thick  is 
sufficient  to  reproduce  the  gravity  low  between  Morenci  and 
Dice  Throw. 

Three  paradigms  were  generated  for  the  Tyrone-Acoma 
line  that  used  the  Morenci-Dice  Throw  thickened  crust,  3.1 
km  and  12.5  km  plutonic  models  as  guides.  Depth  to  major 
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Figure  29.  Gravity  model  along  Jaksha's  (1982)  refraction 
profile,  including  a  3 . 1  km  thick  pluton  (50% 
model) .  Values  for  density  are  given  in  gm/cm^, 
except  for  the  granitic  pluton  which  displays  a 
density  contrast  with  respect  to  its  country 
rock. 
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Figure  30.  Gravity  model  along  Jaksha's  (1982)  refraction 
profile,  including  a  12.5  km  thick  pluton  (20% 
model) .  Values  for  density  are  given  in  gm/cm^, 
except  for  the  granitic  pluton  which  displays  a 
density  contrast  with  respect  to  its  country 
rock. 
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interfaces,  velocity  and  density  were  matched  at  the 
intersection  point  of  these  two  profiles.  That  point  is  in 
Corduroy  Canyon,  approximately  80  km  from  Tyrone  and  125  km 
from  Morenci.  Crustal  structure  and  velocity  were  also 
constrained  near  Tyrone  by  the  previous  work  of  Sinno  and 
others  (1986) .  These  models  were  then  tested  by  gravity  and 
seismic  data  on  the  Tyrone-Acoma  line. 

The  calculated  gravity  results  along  the  Tyrone-Acoma 
line  for  the  two  plutonic  models  are  shown  in  Figures  31  and 
32,  respectively.  The  north  end  of  the  pluton  in  these 
models  is  placed  220  km  from  Tyrone,  the  edge  of  the  sag  in 
gravity  values.  In  these  figures,  the  northern  extent  of 
the  pluton  ,  as  mapped  by  Elston  and  others  (1976a)  and 
suggested  by  Krohn  (1976),  is  150  km  north  of  Tyrone. 

Again,  the  3.1  km  pluton  is  insufficient  to  account  for  the 
entire  gravity  anomaly.  The  12.5  km  pluton  creates  a 
gravity  field  which  is  relatively  close  to  the  observed 
readings.  However,  it  is  clear  from  this  analysis  that  the 
gravity  low  extends  further  north  along  the  profile  than 
predicted  by  Elston  and  others  (1976a) .  Reconciling  this 
problem  requires  extending  the  low-density  body  or 
increasing  the  thickness  of  the  crust  to  the  north. 

The  model  showing  the  anomaly  along  the  Tyrone-Acoma 
refraction  line  due  to  crustal  thickness  alone  is  presented 
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Figure  32.  Gravity  profile  along  the  Tyrone-Acoma  refraction 
line  with  a  12.5  km  thick  buried  pluton.  Values 
for  density  are  in  gm/cro^,  except  for  the 
granitic  pluton  which  displays  a  density  contrast 
with  its  country  rock. 


92 


in  Figure  33.  Based  on  gravity  analysis  alone,  the  crustal 
thickening  and  12.5  km  thick  pluton  models  provide  equally 
acceptable  solutions  to  available  data.  The  validity  of 
each  model  is  tested  by  the  use  of  the  ART  ray-tracing 
program  and  comparing  the  outcomes  to  observed  arrival 
times . 

From  Figure  34,  the  rays  arriving  just  before  the  shift 
in  PmP  arrival  reflect  off  the  Moho  about  100  km  from 
Tyrone.  If  the  travel-time  anomaly  is  caused  by  thinning  of 
the  crust,  it  must  begin  just  beyond  100  km  from  Tyrone.  As 
shown  in  Figure  35,  extreme  thinning  to  23  km  is  required  to 
recreate  this  step.  However,  such  thinning  produces  a 
discrepancy  between  the  predicted  and  observed  gravity 
fields  (Figure  36)  of  almost  50  mgal.  Clearly,  the  gravity 
field  does  not  support  crustal  thinning  of  such  magnitude 
along  this  profile. 

The  time  step  in  the  PmP  arrivals  must  therefore  be 
produced  by  lateral  velocity  variations  within  the  crust.  A 
1.5  sec  shift  in  arrival  times  suggests  a  significant 
feature.  With  an  average  velocity  of  4  km/sec  in  the  near¬ 
surface  layer  and  a  thickness  of  3  km,  the  travel  time  for  a 
vertical  seismic  ray  through  the  section  is  0.75  sec.  If 
this  layer  is  removed  and  replaced  by  upper  crustal  material 
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Figure  33.  Gravity  model  along  the  Tyrone-Acoma  refraction 
line  involving  crustal  thickening  alone.  -lue 
for  density  are  in  gm/cm  . 
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Figure  34.  Seismic  ray-tracing  model  along  the  Tyrone-Acoma 
line  for  a  crust  thickened  to  produce  the 
observed  gravity  anomaly.  Here,  PmP  arrivals  do 
not  recreate  the  observed  time  shift  at  210  km. 
Velocities  are  in  km/sec.  Symbols:  o  =  observed 
arrivals;  +  *  calculated  reflected  arrivals;  x  = 
calculated  refracted  arrivals. 
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Figure  36.  Gravity  profile  along  the  Tyrone-Acoma  refraction 
line  showing  crustal  thinning  to  23  km.  This 
model  is  unsupported  by  the  regional  gravity 
field.  Values  for  density  are  in  gm/cm^. 
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with  a  velocity  of  6  km/sec,  a  difference  in  travel  time  of 
only  0.25  sec  is  expected.  To  create  a  time  shift  of  1.5 
sec,  the  near-surface  layer  must  have  a  velocity  of  1.67 
km/sec,  which  is  unrealistic.  As  a  result,  the  PmP  arrival 
time  step  change  cannot  be  caused  by  variations  in  the  near¬ 
surface  layer  alone,  and  must  be  the  result  of  deeper 
changes  within  the  crust.  This  can  be  accomplished  by 
thickening  the  lower  crustal  layers  in  the  north  at  the 
expense  of  the  upper  crustal  layer,  but  by  doing  so  the 
gravity  field  again  would  be  raised. 

The  presence  of  a  near-surface  plutonic  complex  under 
the  Mogollon  plateau  could  create  the  time  shift  in  PmP.  It 
has  been  mentioned  earlier  that  granitic  batholiths  are  less 
dense  than  their  surroundings.  This  implies  that  seismic 
velocity  within  a  pluton  is  slower  than  the  country  rock. 

The  1.5  sec  step  could  then  indicate  the  point  where  the 
pluton  no  longer  affects  the  ray  path,  suggesting  that  the 
plutonic  complex  extends  approximately  50  km  north  of  Elston 
and  others'  (1976a)  boundary.  This  position  also 
corresponds  to  the  north  end  of  the  sag  superimposed  on  the 
regional  Bouguer  gravity  low  at  that  point. 

The  emplacement  of  such  a  large  plutonic  complex  has 
profound  significance  on  the  overall  structure  of  the  crust. 
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Of  critical  importance  is  the  question:  Where  did  the 
displaced  upper  crustal  material  go? 

Bott  and  Smithson  (1967)  correctly  point  out  that  when 
a  pluton  is  inserted  into  the  upper  crust,  a  like  volume  of 
material  is  displaced  regardless  of  where  the  source  is 
derived.  Forceful,  or  active,  injection  of  fluids  imply 
that  country  rock  is  uplifted,  pushed  aside,  or  is  displaced 
downward  by  stoping.  If  the  primary  displacement  is  upward, 
then  the  20%  model  should  create  a  structural  high  of  12.5 
km  under  the  Mogollon  plateau  relative  to  its  surroundings. 
No  structural  high  of  this  magnitude  has  been  reported  in 
the  region.  If  displacement  is  primarily  outward, 
significant  compressive  features  should  be  seen  in  the  area 
surrounding  the  Mogollon  plateau.  However,  southwestern  New 
Mexico  switched  from  regional  compression  to  tension  during 
the  time  of  emplacement,  and  there  are  no  outward-radiating 
indications  of  stress  from  the  Mogollon  plateau  in  the  mid- 
Tertiary.  Instead,  radially-directed  tension  cracks  occur 
in  a  wide  region  around  the  Mogollon  plateau,  with  the 
center  of  the  pattern  somewhere  to  the  southwest  of 
Magdalena  (Elston  and  others,  1976a) .  Furthermore,  they 
report  that  the  calderas  formed  in  a  neutral  field,  and 
presumably  the  pluton  emplacement  occurred  in  the  same 
field. 
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Therefore,  stoping  must  be  the  dominant  factor  in 
movement  of  country  rock  during  emplacement  of  the  proposed 
mid-Tertiary  pluton  in  southwest  New  Mexico.  Elston  and 
others  (1976a)  propose  that  the  source  region  for  these 
fluids  was  in  the  lower  crust.  Melting  was  due  to  heat  from 
asthenospheric  basalts  that  pooled  at  the  base  of  the  crust. 

With  these  geological  constaints  applied,  the  key 
feature  of  the  final  model  became  apparent:  A  near-surface 
pluton,  extending  through  the  Plains  of  St.  Augustine,  must 
be  present  in  the  crust.  It  must  be  sufficiently  thick  to 
produce  a  fundamental  velocity  change  throughout  the  crust 
via  stoping,  and  account  for  (at  least)  the  sag  within  the 
regional  Bouguer  gravity  low.  It  must  also  satisfy  both  the 
Tyrone-Acoma  and  Jaksha  (1982)  refraction  and  gravity  data. 
The  final  integration  of  these  data  sets  proved  to  be 
difficult  and  time-consuming,  since  perturbing  one  model  had 
a  cascading  effect  on  the  other  three. 

For  the  Tyrone-Acoma  profile,  the  placement  of  the  Moho 
and  upper  mantle  velocities  for  distances  >100  km  from 
Tyrone  were  addressed  by  further  analysis  of  gravity  and 
seismic  data.  Gravity  data  suggested  continued  thickening 
of  the  crust  and/or  pluton.  An  attempt  was  made  to  model  Pn 
data  by  upper  mantle  velocity  changes  alone  as  a  test  for 
pluton  thickening.  Reasonable  velocity  gradients  that 
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yielded  8.0  km/sec  values  at  50  km  were,  however, 
insufficient  to  reproduce  the  Pn  branch  of  the  record 
section  (Figure  37) .  Crustal  thickening  is  therefore 
required  to  account  for  Pn  arrivals  under  the  northern  part 
of  the  pluton.  From  gravity  data,  continued  thickening  to 
the  north  under  the  Colorado  Plateau  is  supported.  This 
implies  that  the  Mogollon-Datil  volcanic  field  currently 
overlaps  the  southeastern  margin  of  the  Colorado  Plateau. 

The  relationship  between  the  magmatism  and  the  southeastern 
margin  during  emplacement  in  mid-Tertiary  time  remains 
unresolved. 

The  bottom  of  the  pluton  in  the  final  model  is 
estimated  to  be  at  10  km,  as  shown  in  Figure  38.  The 
stoping  effect  is  shown  in  this  model  as  a  depression 
recorded  by  the  intermediate  crustal  boundaries.  The 
northern  boundary  of  the  depression  is  more  abrupt  than  that 
to  the  south,  and  may  imply  a  major  fault  of  crustal 
thickness  in  that  region.  However,  no  reflections  from  such 
a  fault  were  observed  in  the  seismic  data,  and  further  work 
is  necessary  to  justify  such  a  conclusion.  A  major  crustal 
fault  and  an  abrupt  sag  under  the  pluton  both  produce  the 
same  effect:  they  bring  faster  velocity  material  nearer  to 
the  surface  north  of  the  boundary.  This  feature,  as  well  as 
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Figure  37.  Ray-tracing  model  for  Tyrone-Acoma  line  ' 
attempting  to  match  Pn  branch  by  upper  mantle 
velocity  anomalies  alone.  Velocities  are  in 
km/sec.  Symbols:  o  =  observed  arrivals;  +  = 
calculated  reflected  arrivals;  x  =  calculated 
refracted  arrivals. 
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Figure  38.  Gravity  profile  for  Tyrone-'Acoma  refraction  line 
final  model.  Values  for  density  are  in  gm/cm^, 
except  for  the  granitic  pluton  which  displays  a 
density  contrast  with  its  country  rock. 
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the  presence  of  the  pluton,  is  necessary  to  create  the 
change  in  arrival  time  of  the  PmP  phase  (Figure  39) . 

To  insure  the  viability  of  the  model  developed  for  the 
Tyrone-Acoma  data,  it  was  applied  to  the  Morenci-Dice  Throw 
profile.  Results  are  shown  in  Figures  40  and  41.  Assuming 
the  east-west  line  is  along  the  strike  of  crustal  structure, 
the  pluton  is  modeled  as  a  slab.  It  is  observed  that  the 
models  are  compatible  in  all  major  respects,  including  the 
reversed  Pn  branch  for  the  Morenci-Dice  Throw  study  (Figure 
42) .  This  model  is  therefore  preferred  over  the  model  of 
crustal  thickening  alone. 

The  model  that  results  from  combined  analysis  of 
seismic  and  gravity  data,  along  with  present  day  knowledge 
of  the  geology  in  the  region,  strongly  supports  a  complex 
history.  The  use  of  both  data  sets  yields  placement  of 
major  crustal  boundaries  with  an  error  estimated  at  +  2  km. 
At  the  north  end  of  the  Tyrone-Acoma  profile,  the  use  of 
gravity  alone  suggests  a  possible  thinning  of  the  crust 
(Figure  38) .  However,  thinning  can  be  avoided  by  increasing 
the  thickness  of  the  upper  and  middle  crust  layers  at  the 
expense  of  the  lower  crust. 

The  time  of  emplacement  of  the  near-surface  plutonic 
complex  is  placed  in  the  mid-Tertiary,  mainly  because  of  the 
assumption  that  it  is  the  source  for  the  overlying 
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Figure  39.  Ray-tracing  model  for  Tyrone-Acoma  line,  final 
configuration.  Velocities  are  in  kin/sec. 
Symbols:  o  =  observed  arrivals;  +  *  calculated 

reflected  arrivals;  x  »  calculated  refracted 
arrivals. 
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Figure  40.  Gravity  profile  for  Jaksha  (1982)  cross-line, 
including  the  final  model.  Values  for  density 
are  in  gm/cm^,  except  for  the  granitic  pluton 
which  displays  a  density  contrast  with  its 
country  rock. 
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Figure  42.  Ray-tracing  model  for  the  Jaksha  (1982)  cross- 
line,  reversed  Pn  branch,  final  configuration. 
Velocities  are  in  km/sec.  Symbols:  o  =  observed 
arrivals;  x  =  calculated  refracted  arrivals. 
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volcanics.  Within  the  constraints  provided  by  these 
geophysical  models,  we  can  conclude  that  the  history  for 
southwestern  New  Mexico  in  mid-Tertiary  time  was  one  of 
large-scale  reorganization  of  the  crust: 

1)  Magmas  formed  from  excess  heat  in  the  lower 
crust  worked  their  way  to  the  near-surface  portion 
of  the  upper  crust.  Movement  must  have  been 
rapid,  as  Elston  and  others  (1976a)  find  little 
evidence  of  crustal  contamination. 

2)  It  is  commonly  held  that  a  basaltic  source  rock 
yields  10%  of  its  volume  to  acidic  fractionation  (Cox 
and  others,  1979) .  The  lower  crust  was  therefore 
either  reduced  in  total  volume,  or  the  loss  of  magma  to 
the  upper  crust  was  replaced  by  cooling  basaltic  fluids 
at  the  base  of  the  crust.  If  the  lower  crust  has  been 
reduced  in  total  volume,  it  is  possible  that  this 
region  was  part  of  the  Colorado  Plateau  prior  to  the 
heating  event.  If  the  fluid  loss  was  replaced  by  new 
magma,  the  Mogollon-Datil  volcanic  field  may  have 
formed  off  the  Colorado  Plateau.  The  solution  to  this 
problem  remains  undetermined  by  this  study. 

3)  The  crustal  column  between  the  near-surface 
pluton  and  the  base  of  the  crust  was  displaced 
downward  by  stoping.  This  implies  that  a  negative 


density  anomaly  was  introduced  at  all  levels  of 
the  crust  in  this  column,  and  slower-velocity 
material  is  pushed  downward  throughout  this 
region. 

4)  Venting  of  the  near-surface  pluton,  especially 
with  caldera  collapse,  provided  a  near-surface 
stoping  effect,  completing  a  total  reorganization 
of  the  crust. 

5)  If  the  initial  heating  event  is  caused  by  the 
introduction  of  asthenospheric  fluids  at  the  base 
of  the  crust,  it  is  implied  that  the  upper  mantle 
section  of  the  lithosphere  has  been  disturbed  in 
addition  to  the  crust.  Presumably  these  magmas 
came  from  depth,  and  implies  a  stoping  effect  in 
the  upper  mantle.  Then  the  mid-Tertiary  events  in 
southeastern  Arizona  and  southwestern  New  Mexico 
actually  represent  reorganization  of  the  entire 
lithosphere. 


THERMAL  MODELING  OF  THE  LITHOSPHERE, 

ARIZONA  AND  NEW  MEXICO 

The  present  structure  of  the  crust,  along  with  known 
geology,  suggests  a  complex  tectonic  history  in  southeastern 
Arizona  and  southwestern  New  Mexico.  Although  the  nature 
and  chemistry  of  the  volcanism  has  changed  over  time,  the 
continued  activity  indicates  heating  at  depth  at  least  since 
the  mid-Tertiary.  The  model  describing  present  day  crustal 
structure  as  derived  in  the  above  section  must  be  understood 
with  respect  to  preceding  events. 

The  complex  nature  of  observed  heat  flow  in  the 
Colorado  Plateau,  Basin  and  Range  and  Rio  Grande  Rift 
provinces  and  their  transition  zones  is  introduced  in  a 
previous  section  in  this  report.  However,  a  brief  review  is 
necessary  for  the  present  discussion.  This  section 
describes  a  two-dimensional  temperature  model  that  attempts 
to  recreate  the  present-day  lithospheric  thetToal  structure 
by  a  finite  difference  method.  In  addition,  it  calculates 
heat  flow  and  isostatic  uplift  over  the  model.  By 
comparison  to  modern  heat  flow  and  elevation  the  input 
parameters  can  be  constrained,  and  reasonable  physical 
paradigms  can  be  investigated. 
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Heat  flow  values  for  Arizona  and  New  Mexico  were 
extracted  from  a  national  data  base  (P.  Morgan  and  J.  Sass, 
personal  communication,  1988)  and  are  plotted  in  Figure  20. 
This  map  demonstrates  that  high  heat  flow  values  are 
associated  with  the  boundary  of  the  Colorado  Plateau.  This 
relationship  implies  that:  1)  a  thermal  pulse  is  spatially 
associated  with  the  boundary,  i.e.  in  the  Basin  and  Range 
and  Rio  Grande  Rift  (Seager  and  Morgan,  1979) ;  or  2)  faults 
or  other  conduits  allowing  hot  fluids  to  be  transported  to 
the  surface  are  associated  with  the  plateau  margin 
(Swanberg,  1979) . 

The  Colorado  Plateau  displays  heat  flow  values  near  or 
slightly  above  those  found  in  stable  continental  regions 
(Blackwell,  1978;  Lachenbruch  and  Sass,  1978;  Bodell  and 
Chapman,  1982;  Reiter  and  Clarkson,  1983).  This  heat  flow 
is  highly  variable,  however,  and  care  must  be  taken  in  the 
interpretation  of  available  data  (Reiter  and  Clarkson, 

1983) .  Shearer  and  Reiter  (1981)  and  Sass  and  others  (1982) 
discuss  the  influence  of  groundwater  circulation  on  heat 
flow  values  in  the  San  Francisco  volcanic  field  on  the 
plateau  margin  in  Arizona.  Significant  hydrothermal  effects 
are  also  noted  in  the  Rio  Grande  Rift  (Reiter  and  others, 
1978)  .  This  concern  is  typical  of  the  problems  in  much  of 
the  western  United  States  (Blackwell,  1978;  Sass  and  others, 
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1982)  .  Figure  43  shows  that  the  transition  between  the 
Colorado  Plateau  and  the  extensional  provinces  to  the  south 
and  east  is  complex.  It  appears  to  be  at  least  100  km  wide 
in  Arizona  and  New  Mexico  (Bodell  and  Chapman,  1982; 

Swanberg  and  Morgan,  1985) . 

The  recent  uplift  of  the  region  indicates  some  heat  has 
been  taken  up  by  thermal  expansion.  However,  Morgan  and 
Swanberg  (1985)  and  Hinojosa  (1989)  show  that  thermal 
expansion  is  insufficient  for  the  total  amount  of  uplift 
under  the  Colorado  Plateau.  In  fact,  Morgan  (1983)  states 
that  a  surface  heat  flow  anomaly  should  not  develop  until 
after  lithospheric  thinning  has  stopped.  It  has  been  shown 
in  a  previous  section  that  geophysical  data  indicate  that 
the  upper  mantle  is  hot  under  the  plateau,  even  though  the 
surface  heat  flow  indicates  lower  temperature  at  depth. 

This  tends  to  support  Morgan's  (1983)  study,  but  does  not 
address  the  present  elevation  and  uplift  of  the  plateau. 

The  causes  and  consequences  of  the  geophysical  data  and 
present  day  tectonic  state  of  the  Colorado  Plateau  are 
beyond  the  scope  of  this  study. 

Heat  flow  values  have  been  separated  into  "deep"  values 
and  "shallow"  values  within  the  Basin  and  Range  province  in 
Arizona  (Shearer  and  Reiter,  1981)  and  the  Colorado  Plateau 
(Reiter  and  Clarkson,  1983).  In  Arizona,  temperature 
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gradients  measured  at  depths  >650  m  ("deep")  yield  more 
consistent  values  (i.e.  smaller  variance)  than  the  province 
as  a  whole.  The  heat  flow/heat  generation  regression  for 
deep  data  yields  a  reduced  heat  flow  (intercept)  of  1.71  HFU 
for  the  Arizona  Basin  and  Range.  The  "shallow"  heat  flow, 
determined  from  well  data,  is  1.94  +  0.12  HFU.  The  low 
standard  error  is  suggested  to  indicate  a  uniform  radiogenic 
layer,  which  in  turn  implies  a  homogeneous  crust (Shearer  and 
Reiter  1981) .  For  the  Colorado  Plateau,  data  from  depths 
>750  m  also  display  much  lower  lateral  variability  than 
near-surface  measurements  (Reiter  and  Clarkson,  1983) . 
Furthermore,  Reiter  and  Clarkson  find  a  higher  heat  flow 
mean  in  the  deep  data  (1.53  HFU)  than  in  shallow 
measurements  (1.31  HFU).  Greater  lateral  consistency  in  the 
deep  values  imply  an  equal  radiogenic  contribution  to  the 
heat  flow,  and  also  consistent  upper  mantle  temperatures. 
Within  the  limits  of  statistical  analysis,  Shearer  and 
Reiter  (1981)  show  that  the  difference  between  the  mean  heat 
flow  of  the  Colorado  Plateau  and  Basin  and  Range  provinces 
may  be  as  little  as  0.3  to  0.5  HFU.  This  implies  that  their 
tectonic  history  is  more  similar  than  previously  thought, 
and  that  the  crust  has  responded  differently  in  each 
province.  However,  the  deep  data  are  sparse  in  both 
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provinces,  and  conclusions  based  on  them  should  be  viewed 
with  caution. 

The  Rio  Grande  Rift  and  Basin  and  Range  provinces 
display  high  heat  flow  near  faults  and  volcanic  centers,  and 
intermediate  to  slightly  elevated  heat  flow  elsewhere 
(Reiter  and  others,  1975,  1978;  Shearer  and  Reiter,  1981; 
Reiter  and  Clarkson,  1983).  However,  Swanberg  (1979)  states 
that  the  Rio  Grande  Rift  displays  consistently  high  heat 
flow  values  along  its  axis  (Figure  44) .  High  heat  flow  in 
the  rift  is  thought  to  be  caused  by  tectonic  and  magmatic 
sources  in  the  crust  and  upper  mantle  (Edwards  and  others, 
1978).  Heat  flow  values  exceeding  3.0  HFU  appear  to  the 
west  of  El  Paso,  but  are  not  found  near  the  Colorado  Plateau 
in  New  Mexico. 

For  a  continental  rift  with  dimensions  near  those  of 
the  East  African  Rift  in  Kenya,  surprisingly  little 
volcanism  has  occurred  within  the  Rio  Grande  structure 
(Baldridge  and  others,  1984).  Magma  compositions  in  and 
near  the  rift  are  highly  variable  and  do  not  necessarily 
reflect  their  tectonic  setting.  Sources  are  also  varied  in 
depth  of  origin  for  the  extrusives  found  in  and  around  the 
rift  (Baldridge  and  others,  1987;  Perry  and  others,  1987, 
1988;  Elston,  1989) .  Because  of  the  small  total  volcanic 
volume,  however,  Baldridge  and  others  (1984)  conclude  that 
no  major  thermal  pulse  is  present  in  the  upper  mantle 
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Figure  44. 


Heat  flow  along  the  Rio  Grande  Rift  and  adjacent 
areas.  After  Swanberg  (1979). 
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beneath  the  Rio  Grande  Rift.  Although  pooling  at  mid- 
crustal  depths  appears  to  occur  along  the  rift  (e.g.  Sanford 
and  others,  1977) ,  the  volume  of  such  bodies  appears  to  be 
minor  or  is  undetectable  by  available  heat  flow  data. 

Although  the  highest  heat  flow  values  in  the  Rio  Grande 
Rift  occur  near  deep  faults  and  volcanic  centers,  the  higher 
heat  flow  trend  along  its  axis  (Swanberg,  1979)  indicates  a 
more  fundamental  thermal  disturbance.  The  lack  of  extensive 
volcanism  and  plutonism  (Baldridge  and  others,  1984)  may  be 
the  result  of  a  crustal  composition  that  resists  melting  at 
this  time.  Sparse  earthquake  activity  (Jaksha  and  Sanford, 
1986)  implies  little  voluminous  movement  of  magma  in  the 
crust.  For  this  study,  it  is  assumed  that  a  thermal 
disturbance  is  present  in  the  upper  mantle,  and  that  the 
small  volcanic  volume  indicates  the  primary  transport  of 
thermal  energy  is  conductive. 

A  two-dimensional  finite  difference  conductive  heat 
transfer  program  has  been  created  in  an  attempt  to  analyze 
the  phenomena  described  above.  Its  derivation  is  presented 
in  the  following  section,  and  the  resulting  program  HF2D  is 
presented  in  Appendix  V.  It  is  applied  as  a  cross-section 
of  the  Colorado  Plateau  boundary,  and  is  assumed  to  be 
infinite  in  the  y-direction.  The  origin  of  the  model  is 
placed  at  the  physical  axis  of  high  heat  flow  readings  (z  = 
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0) ,  and  temperature  perterbations  are  symmetric  about  the  z- 
axis  of  the  model  in  the  x-direction. 
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Derivation  of  the  Algorithm 

The  analytical  form  of  the  two-dimensional  heat  flow 
equation  (HF2D)  is: 

11  =  K  /si  +  K'i 

6t  V  5x*  5 1  > 

Where  T  =  Temperature 
t  =  Time 

K  =  Thermal  diffusivity 
X  =  Distance  in  horizontal  direction 
2  =  Distance  in  vertical  direction 


It  is  desirable  that  the  discrete  form  of  this  equation 
be  unconditionally  stable.  An  equation  that  satisfies  this 
form  is  the  alternating-direction  implicit  (ADI)  method 
described  by  Anderson  and  others  (1984) .  This  equation 
splits  the  forward  step  in  time  into  two  "half-steps", 
taking  the  following  form: 


Step  2 


At/2 
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Where 


At 

6  I 


Ti  ♦  i,  i 


■  • 

2T|,j  +  Ti_i,4 


(  AX)* 


a.  “ 

5*  T, 


» .  J  - 


(  Az)* 


n  =*  Time  from  which  this  step  is  being  made 

n+1/2  -  Intermediate,  or  "half-step” 
n+1  =  Next  time  step 

i,j  *  X  and  z  points  on  a  grid  of  dimension 
N  X  N 

This  method  is  second  order  accurate  with  a  truncation 
error  of:  0  (  (At)*,  (Ax)*  and  (Az)*  1  • 


for  step  1, 

Ti.j  -  T,,,  =  KAt  -  2T,,,  ♦  T,.i.,  '  -2T,,,  + 

-( - + - ^ - 

2  (Ax)*  (Az) 


=  K  At 

2(  Az)2 


Let 


a  = 


K  At _ 

2  ( A  X)^ 


and 


b 
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Then  we  find  that  for  step  1, 


■*l/t  a*!/!  m*l/2 

+  (l-f2a)T|,j  -aT|_i,4  ~ 


+  (l-2b)Ti,j  +  bTj.j.i 


For  step  2 


i*l/t  a*!/! 


Tj.^  -  T, ■;;■*=  a(TMi.i  -2Tr;;'v  t’:;:')  +  b(T,:j:,  -2Ti'V  ♦  t.V/.j) 


ft* i  ••  I 

■bT,,jM  +  (l  +  2b)T,,4 


bT,:;!, 


aTu,‘.7  +  (i-2a)Tr;;'*  +  aT,*:;;; 


For  the  present  discussion,  we  will  apply  the  following 
boundary  conditions: 


fix  X  «  0,L 


Heat  flow  through  the  lateral 
edges  is  set  =  0 


Temperature  at  the  surface  is 


z  =  0 


set  =  0 


z  «  H 


=  f(x)  Temperature  at  the  base  is 

some  function  of  x 


I 
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Numerically,  the  lateral  change  in  temperature  at  any  point 
can  be  defined  as: 

5  X  2  A  X 

since  the  heat  flow  =  0,  then  j  =  T,.,  ^  .  This  step  is 


necessary  for  step  1,  when  the  points  Tj  j  and  j  are 
evaluated.  The  matrix  notation  for  the  first  step  is: 

n+1/2 


I 

I 
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Where  the  superscripts  n  and  n+1/2  refer  to  the  time  step 
for  these  temperature  values.  This  equation  increments  i 
from  1  to  N  while  holding  j  constant.  After  the  i-loop 
terminates,  j  increases  by  1  and  the  process  repeats.  Since 
Tj  ^  and  are  held  fixed  (upper  and  lower  boundary 
conditions) ,  they  are  not  necessary  in  the  outer  loop.  Thus 
the  matrix  equation  takes  the  form  A  ^  =  y  where  A  is  an 
NxN  matrix,  x  is  the  Nxl  vector  of  temperatures  for  row  j  at 
time  step  n+1/2,  and  “y  is  the  Nxl  vector  of  temperature- 
spatial  parameters  at  time  step  n.  Noting  that  A  is  a 
tridiagonal  matrix,  much  computer  time  can  be  saved  in 
calculating  A*’  by  using  program  TRIDAG  (Press  and  others, 
1986;  also  see  Appendix  V). 

The  second  step  requires  that  i  be  held  constant  while 
j  increments  form  2  to  N-1  (For  the  same  reason  as  in  step 
1) .  However,  the  right-hand  side  of  equation  (3)  depends  on 
the  value  of  i.  For  i  =  1  and  i  =  N,  the  heat  flow  =  0 
boundazry  condition  must  be  applied  for  Tq  j  and  T^j^,  j . 
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Thus: 

For  i  =  1 


n+r 


Note  that  the  left-hand  matrix  is  no  longer  square,  it  is  of 
dimension  (N-2)xN.  Also,  the  right-hand  side  combines  to 
form  an  (N-2)xl  vector.  However,  it  is  noted  that  the  top 
and  bottom  elements  in  the  t”*’  vector  (left-hand  side)  are 
points  at  which  the  temperature  never  changes  (boundary 
conditions) .  Thus  the  a  matrix  can  be  made  square  and 
tridiagonal,  and  the  right-hand  side  becomes  Nxl,  by  the 
following  modifications: 
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1  0  0  .  .  •  0 

Ti,r 

-b  ( l+2b)  -b  .  • 

Ti,t 

0  -b  {l+2b)  -b  .  •  • 

Ti,. 

.  a  e  .  •  •  • 

• 

Ti.4 

.  a  •  •  *  *  9 

• 

-b  (l+2b)  -b 

• 

0  a  a  .  0  0  1 

Ti,« 

n+1/2 

Tui 

2aT,,,  ♦  (l-2a)Tx,, 

2aT,.,  ♦  (l«2a)Ti,, 

2aT,.4  ♦  (l-2a)T».« 

• 

2aTf,ii_i  +  ( l*”2a )  Ti,  n-i 

T,.« 


Again,  this  is  in  the  form  5  7-7  where  A  is  NxN,  and  x 
and  7  are  both  Nxl.  A  similar  derivation  for  i  =  N  yields: 
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n+1 
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• 

-  • 

0 

0 
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Tii,  i 

(l-2a)T«,,  +  2aT,-i., 
(l-2a)T|i,j  ♦  2aT||.t,s 
(l-2a)T*,4  +  2aT,.i.4 


n+l/2 


(l-2a)T«,,.|  +  2aT*.x.,-i 

T«.. 
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Finally,  for  i  from  2  to  N~l: 


1  0  0  .  •  •  ° 

—  “1 
Ti.i 

-b  {l+2b)  -b  .  •  •  * 

Ti,. 

0  -b  (l+2b)  -b  .  •  • 

T,., 

. 

• 

T».4 

. 0 
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o 

• 
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+ 
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2aT2, 1 

+ 
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2aT2,4 
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Since  i  is  held  constant  while  j  increments,  this  calculates 
one  row  at  a  time.  The  result  is  the  2D  temperature  profile 
at  time  n+1. 


Consequences  of  Lithospheric  Heating 

An  increase  in  temperature  at  the  base  of  the 
lithosphere  appears  at  the  surface  after  a  certain  lag  time, 
which  is  a  function  of  the  average  thermal  diffusivity  of 
the  rock  column.  If  the  transport  of  heat  is  conductive, 
the  result  will  be  an  increase  in  surface  heat  flow.  In 
addition,  an  increase  in  temperature  will  result  in  a 
thermal  expansion  of  the  rock  material.  This  expansion 
changes  the  density  of  the  material,  with  the  change  being  a 
function  of  the  coefficient  of  thermal  expansion.  As  a 
result,  an  increase  in  elevation  in  response  to  isostatic 
uplift  can  be  calculated.  Program  HF2D  has  been  designed  to 
yield  these  quantities  for  a  homogeneous  two-dimensional 
model . 


Surface  Heat  Flow 

The  surface  heat  flow  is  the  product  of  the  geothermal 
gradient  and  the  thermal  conductivity  of  the  rock  material. 
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To  obtain  the  gradient  at  a  given  horizontal  position,  the 
program  takes  the  difference  between  the  top  two  points  in 
the  model  divided  by  the  distance  between  them.  The  result 
is  a  set  of  heat  flow  values  as  a  function  of  distance  from 
the  z-axis  in  the  model. 

Isostatic  Uplift 

The  isostatic  uplift  at  a  given  point  is  the  integrated 
change  in  density  along  the  underlying  rock  column.  The 
integration  algorithm  chosen  for  use  in  this  program  is  the 
1/3  Simpson's  rule,  determined  for  an  odd  number  of  points. 
The  result  is  a  set  of  values  for  uplift  in  kilometers  as  a 
function  of  distance  from  the  z-axis  in  the  model. 

Consequences  of  the  Model 

Within  certain  restrictions  discussed  below,  any  of  the 
data  sets  calculated  in  this  program  can  be  compared  to  the 
present  day  values  found  in  New  Mexico  and  Arizona.  In  this 
way,  each  provides  constraints  on  the  amount,  extent  and 
duration  required  for  a  thermal  pulse  at  the  base  of  the 
lithosphere.  However,  these  values  taken  together  and 
compared  to  present  day  observations  of  heat  flow  and 
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elevation  will  provide  an  optimum  constraint  to  determine 
the  Cenozoic  thermal  history  of  Arizona  and  New  Mexico. 

Discussion 

A  series  of  thermal  perturbations  have  been  applied  to 
a  steady-state,  stable  continental  temperature  gradient  in 
order  to  determine  a  reasonable  thermal  history  in  the  area. 
The  present  day  heat  flow  anomaly  associated  with  the  Rio 
Grande  Rift  can  be  viewed  as  symmetric  over  a  significant 
distance  along  its  axis  in  New  Mexico  (Figure  44) .  Because 
the  surface  heat  flow  falls  off  towards  both  the  Colorado 
Plateau  and  Great  Plains  provinces,  the  perturbation  is  in 
the  form  of  a  gaussian  'bell'  curve  described  by: 

T  =  C  X  exp  (-a  X  d^) 

Where  T  =  Temperature 

C  =  Amplitude  of  curve  at  distance  x  =  0 
a  =  Power  coefficient,  which  determines  the  width 
of  the  bell 

d  B  Offset  distance  in  km. 

In  order  to  minimize  spatial  error  in  this  analysis, 
only  half  of  the  bell  curve  is  used  in  the  calculation.  The 
x-values  can  then  be  placed  twice  as  close  as  if  the  whole 
bell  curve  is  used.  This  is  allowable  because  of  the 
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symmetry  inherent  in  the  problem.  The  symmetry  also  implies 
that  the  lateral  heat  flow  across  the  axis  of  the  bell  is 
zero,  a  condition  required  for  the  mathematical  derivation. 
For  the  present  analysis,  therefore,  the  axis  of  the 
gaussian  perturbation  coincides  with  the  axis  of  high  heat 
flow  readings  along  the  Rio  Grande  Rift  in  central  New 
Mexico.  The  program  assumes  that  the  thermal  pulse  is 
infinite  in  extent  along  the  line  perpendicular  to  the  bell 
curve. 

These  thermal  pulses  have  been  applied  at  both  100  km 
and  30  km  depth  to  calculate  their  effect  on  heat  flow  and 
isostatic  uplift.  However,  crustal  thinning  has  taken  place 
both  in  the  Basin  and  Range  province  (Gish  and  others,  1981; 
Sinno  and  others,  1981;  this  study)  and  Rio  Grande  Rift 
(Daggett  and  others,  1986;  Sinno  and  others,  1986;  this 
study) .  Since  the  lower  crust  has  been  replaced  by  mantle 
(asthenosphere  ?)  material,  a  convective  component  has  been 
added  in  those  regions,  which  is  not  addressed  by  program 
HF2D.  This  has  a  direct  effect  on  the  isostatic  uplift 
calculations,  which  are  based  on  thermal  expansion  alone. 
Therefore,  for  the  cases  in  which  heating  is  confined  to  the 
base  of  the  lithosphere  (100  km  depth) ,  the  calculated 
values  for  isostatic  uplift  are  not  correlatable  to 
observations  from  the  Basin  and  Range  and  Rio  Grande  Rift  in 


132 


Arizona  and  New  Mexico.  However,  since  the  entire  region 
has  been  uplifted,  the  isostatic  uplift  calculated  by  HF2D 
may  have  a  broad  meaning  to  the  discussion.  That  is,  uplift 
has  occurred  over  the  entire  region,  affecting  all  three 
provinces.  The  calculation  of  isostatic  uplift  can 
therefore  be  used  as  a  general  constraint  to  subsurface 
temperatures  and  their  effect  on  surface  heat  flow.  Also, 
all  of  the  values  created  in  the  program  should  have 
quantifiable  meaning  for  perturbations  at  30  km,  since  they 
involve  only  the  crust. 

The  initial  estimates  for  temperature  and  density  are 
similar  for  all  models  run  in  this  study,  and  are  presented 
in  Figures  45  and  46  respectively.  The  initial  steady-state 
temperature  gradient  includes  a  radioactive  contribution  for 
the  top  10  km,  underlain  by  a  linear  increase  in  temperature 
to  the  bottom  of  the  lithosphere.  After  the  program  started 
calculating  temperature  evolution,  the  near  surface 
radioactive  contribution  was  removed.  This  resulted  in  an 
effective  downwarp  in  elevation  and  lowering  of  the  surface 
heat  flow  over  the  entire  model,  which  was  overcome  by 
increasing  the  thermal  pulse  at  depth.  In  this  way,  the 
maximum  thermal  perturbation  could  be  calculated  for  the 
observed  uplift  and  heat  flow  values.  The  density 
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Figure  45.  Initial  temperature  distribution  for  program 

HF2D.  This  is  a  typical  steady-state  temperature 
gradient  with  a  10  km  thick  radioactive  layer  in 
the  upper  crust,  underlain  by  nonradioactive 
material . 
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Figure  46.  Initial  density  distribution  for  program  HF2D. 

The  upper  40  km  consist  of  typical  crustal 
density  values,  and  the  lower  60  km  is  that  of 
mantle  material. 
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distribution  is  for  a  crust  of  40  km  thickness,  with  a 
gradient  resulting  from  overburden  pressure. 

The  amplitude  and  width  of  the  temperature  pulse 
directly  affects  the  resulting  heat  flow,  uplift  and  gravity 
readings.  Two  curves  of  different  width  were  used 
throughout  this  analysis:  1)  using  a  power  exponent  of 
0.0002,  a  curve  with  half  width  =  60  km  was  generated  and  is 
displayed  with  a  zero  offset  amplitude  of  100°C  in  Figure 
47;  and  2)  using  the  same  amplitude,  Figure  48  shows  a  half¬ 
width  of  37.5  km  that  reflects  a  power  coefficient  of 
0.0005.  By  varying  the  amplitude,  these  perturbations  are 
allowed  to  progress  into  the  model  with  time,  and  their 
results  are  evaluated  with  respect  to  observed  values. 

A  cross-section  of  the  heat  flow  in  the  Rio  Grande  Rift 
is  provided  by  Reiter  and  others  (1978) .  They  show  an 
average  peak  value  of  2.7  HFU  over  the  axis  of  the  anomaly, 
falling  off  to  approximately  1.8  HFU  about  90  km  away 
(Figure  49) .  If  the  thermal  anomaly  lies  at  100  km  depth, 
the  heat  flow  values  generated  by  HF2D  approach  their  curve 
at  the  axis  only  by  a  1500°C  amplitude  with  a  power 
coefficient  of  0.0002  (Figure  50).  The  "background"  heat 
flow  observed  at  distance  is  1.5  HFU,  which  is  less  than 
Reiter  and  others  (1978)  background  of  1.8  HFU.  The 
difference  is  caused  by  dissimilar  steady-state  temperature 
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Figure  47.  Gaussian  curve  with  amplitude  =  100,  exponential 
coefficient  term  ®=  0.0005.  Half-width  of  this 
pulse  =  75  km  due  to  symmetry  about  the  origin. 
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Figure  48.  Gaussian  curve  with  amplitude  =  100,  exponential 
coefficient  term  *  0.0002.  Half-width  of  this 
pulse  *  120  km  due  to  symmetry  about  the  origin. 
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Figure  49. 


Thermal  model  of  Reiter  and  others  (1978)  showing 
heat  flow  rise  over  the  Rio  Grande  Rift. 
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Figure  50.  Calculated  heat  flow  for  a  thermal  pulse  at  100 
km  depth  after  40  Ma.  The  pulse  is  a  gaussian 
curve  centered  on  the  origin,  with  a  1500  C 
amplitude  and  power  coefficient  of  0.0002.  The 
heat  flow  generated  by  this  model  is  similar  to 
that  found  in  Reiter  and  others  (1978).  As  in  all 
following  figures,  the  origin  is  coincident  with 
the  axis  of  heat  flow  along  the  Rio  Grande  Rift. 
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gradients.  It  is  observed  that  the  heat  flow  anomaly  is  too 
broad  for  the  rift.  Furthermore,  the  calculated  uplift  of 
5.5  km  on  the  z-axis  (Figure  51)  shows  this  model  to  be 
unreasonable,  since  the  regional  uplift  of  the  Colorado 
Plateau  is  observed  to  be  about  2  km  (Morgan  and  Swanberg, 
1985)  .  A  perturbation  of  1500°C  at  100  km  depth  produces 
temperatures  in  excess  of  2700°C,  yielding  a  temperature 
gradient  of  27°C/km  through  the  entire  lithosphere  on  the 
axis  for  the  steady-state  result.  This  value  is 
unacceptably  high  for  temperatures  ai-  100  km  depth. 
Furthermore,  steady-state  temperatures  had  not  been  achieved 
even  after  40  Ma,  equivalent  to  the  time  from  the  Laramide 
to  the  present.  For  these  reasons,  this  model  has  been 
rejected  as  a  solution  to  the  present  day  tectonic  state  of 
the  lithosphere. 

Further  analysis  of  Figure  51  determines  that  away  from 
the  axis,  the  values  for  isostatic  uplift  become  negative 
and  indicate  a  downwarp  of  about  350  m.  This  is 
attributable  to  the  loss  of  the  radioactive  component  to 
temperature  values  during  the  model  evolution.  Since  the 
thermal  effects  of  a  temperature  pulse  at  the  axis  should 
have  little  effect  at  distance,  the  effect  of  removing  the 
radioactive  component  is  that  of  a  shift  of  the  entire 
profile  -350  m.  The  correct  uplift  is  therefore  estimated 
to  be  nearly  6  km  at  the  axis  for  the  model  presented  above. 


Isostotic  Uplift  km 


Figure  51.  Calculated  isostatic  uplift  for  perturbation 

described  in  Figure  50.  The  negative  values  at 
large  distances  are  caused  by  loss  of  the  radio¬ 
active  component  during  temperature  evolution, 
which  results  in  a  DC  shift  of  -  0.35  km  to  the 
data.  The  correct  uplift  is  nearly  6  km. 
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The  correct  values  for  similar  figures  are  hereafter 
referred  to  as  uplift  after  "radioactive  correction". 

It  is  reasonable  to  expect  the  surface  heat  flow  to 
decline  from  loss  of  the  radioactive  component  as  the  model 
evolves.  However,  with  an  initial  near-surface  temperature 
gradient  of  24°C  and  conductivity  of  2.5  W/m°C,  the  initial 
heat  flow  is  1.5  HFU.  The  lack  of  change  of  heat  flow  over 
time  indicates  that  isostatic  uplift  is  more  sensitive  to 
temperature  changes  throughout  the  lithosphere.  The  sole 
use  of  heat  flow  values  can  therefore  be  misleading  when 
estimating  tectonic  change  in  the  lithosphere. 

Further  evidence  for  this  problem  is  presented  in 
Figures  52  and  53,  in  which  the  thermal  perturbation  at  100 
km  depth  has  been  reduced  to  400‘’c.  The  heat  flow  (Figure 
52)  falls  to  1.75  HFU  on  the  axis  representing  a  33%  drop  in 
value.  The  total  uplift  after  the  radioactive  correction 
has  been  reduced  to  approximately  1.5  km  (Figure  53),  a  drop 
of  75%.  In  this  case,  both  the  heat  flow  and  the  uplift 
values  are  too  low  when  compared  to  observations  in  central 
and  southern  New  Mexico.  For  this  example,  the  surface  heat 
flow  observed  near  the  border  of  the  Colorado  Plateau  cannot 
be  created  by  a  single  thermal  pulse  at  the  base  of  the 
lithosphere. 


Heat  Flow,  HFU 


Distance.  I<m 


Figure  52.  Calculated  heat  flow  for  a  400°C  thermal  pulse  at 
100  Ian  depth.  Power  coefficient  =  0.0002. 


Distance,  km 


Figure  53.  Calculated  isostatic  uplift  for  a  400°C  thermal 
pulse  at  100  km  depth.  Power  coefficient  = 
0.0002.  Total  uplift  after  radioactive 
correction  is  1.5  km  (see  Figure  51). 
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The  effects  of  a  thermal  pulse  at  the  base  of  the  crust 
are  therefore  investigated.  Because  the  highest  heat  flow 
readings  are  adjacent  to  the  Colorado  Plateau  but  within  the 
Basin  and  Range  and  Rio  Grande  Rift  provinces,  the  base  of 
the  crust  is  placed  at  30  km.  A  350°C  pulse  on  the  axis 
with  a  power  coefficient  =  0.0005  yields  a  heat  flow  of  2.8 
HFU  above  the  center  of  the  anomaly  (Figure  54) ,  which  is 
close  to  observed  readings  along  the  rift.  The  heat  flow 
bulge  also  shows  a  half-width  of  about  40  km,  which  agrees 
reasonably  well  with  the  data  of  Reiter  and  others  (1978) . 
However,  this  thermal  pulse  yields  an  uplift  after 
radioactive  correction  of  only  0.39  km  at  the  axis  (Figure 
55)  indicating  that  this  pulse  is  insufficient  for  the 
observed  regional  uplift. 

At  steady-state,  this  pulse  would  create  a  25°C/km 
gradient  to  the  anomaly  at  the  base  of  the  crust.  From  the 
arguments  in  the  previous  examples,  this  gradient  cannot 
extend  to  the  base  of  the  lithosphere.  The  gradient  must 
therefore  flatten  below  the  crust  in  order  to  prevent 
excessive  uplift.  Steady-state  temperatures  are  reached  for 
this  perturbation  in  less  than  5.5  Ma,  which  is  close  to  the 
time  of  the  most  recent  Basin  and  Range  type  extension  and 
volcanism  in  the  area,  including  the  Jemez  lineament  (Luedke 
and  Smith,  1978) . 


Heal  Flow.  HFU 


Distance,  km 


Figure  54.  Calculated  heat  flow  for  a  350°C  thermal  pulse  at 
30  km  depth. 


Distance,  km 


Figure  55.  Calculated  isostatic  uplift  for  a  350°c  thermal 
pulse  at  30  km  depth.  Total  amount  of  uplift 
after  radioactive  correction  is  0.39  km  (see 
Figure  51) . 
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The  cause  of  this  temperature  pulse  may  be  the 
introduction  of  basaltic  magma  at  the  base  of  the  crust,  as 
proposed  by  Elston  and  others  (1976a)  and  Davis  and  others 
(1989) .  From  that  position,  the  magma  may  have  become  the 
source  for  recent  basaltic  volcanism  in  the  region  (Luedke 
and  Smith,  1978) .  Its  presence  may  also  have  caused  the 
recent  thermal  erosion  of  the  crust  described  by  Sinno  and 
others  (1981) . 

Because  of  the  discrepancies  between  isostatic  uplift 
and  heat  flow  for  models  based  on  temperature  pulses  at  100 
km  and  30  km,  they  are  insufficient  to  explain  present-day 
values  and  therefore  tectonic  history.  A  better  solution 
may  be  found  by  a  combination  of  the  above  models.  It  is 
clear  from  these  examples  that  temperature  anomalies  at  lOO 
km  depth  have  greater  influence  on  uplift,  and  perturbations 
at  30  km  have  more  impact  on  heat  flow.  By  superimposing 
the  400°C  temperature  at  100  km  and  the  350°C  pulse  at  30 
km,  a  solution  that  satisfies  observed  data  may  be  obtained. 

From  arguments  presented  above,  the  temperature 
gradient  must  flatten  below  the  crust  if  a  350°C  pulse  is 
located  at  30  ]cra  depth.  Since  heat  flow  is  dependent  on 
temperature  gradient,  the  deep  anomaly  has  no  effect  for 
this  condition.  Therefore,  no  difficulties  are  encountered 
for  heat  flow  data  in  this  model. 
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Isostatic  uplift  must  be  some  combination  of  the  two 
models.  In  effect,  raising  the  temperature  by  400°C  at  100 
km  depth  raises  the  thermal  gradient  by  4°C/kro,  At  3  0  km, 
the  temperature  is  raised  120°C.  The  uplift  associated  with 
this  pulse  after  radioactive  correction  is  about  1.5  km. 
Raising  the  temperature  350°C  at  30  km  depth  with  respect  to 
initial  conditions  corresponds  to  an  actual  increase  of  only 
230°C  over  the  new  "background"  temperature.  The  additional 
uplift  for  the  crustal  component  after  radioactive 
correction  is  about  0.23  km  on  the  z-axis  (Figure  56) ,  and 
the  extra  uplift  for  increased  temperatures  from  30  km  to 
100  km  is  about  0.64  km.  The  total  uplift  is  about  2.37  km, 
which  is  slightly  higher  than  the  minimum  amount  of  uplift 
required  by  Morgan  and  Swanberg  (1985) .  With  heat  flow 
data,  the  combination  of  these  temperatures  at  30  km  and  100 
km  depth  appears  to  satisfy  observed  data.  The  steady-stat® 
temperature  gradient  is  achieved  for  the  deep  pulse  at  27.8 
Ma  implying  that  it  was  initiated  during  the  mid-Tertiary. 

The  values  for  heat  flow  and  isostatic  uplift 
calculated  from  program  HF2D  are  based  on  some  assumptions 
that  are  violated  by  this  explanation.  For  example,  pooling 
of  basalt  fluids  at  the  base  of  the  crust  implies  a 
significant  convective  component  to  heat  flow,  while  HF2D 
assumes  heat  transfer  by  conduction  alone.  Furthermore, 


Distance,  km 


Figure  56.  Calculated  isostatic  uplift  for  a  230°C  thermal 
pulse  at  30  km  depth.  Total  uplift  after  radio 
active  correction  (see  Figure  51)  is  0.25  km. 
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thinning  of  the  crust  introduces  mantle  material  as  a 
replacement  for  crustal  material,  which  violates  the 
program's  assumption  of  change  in  density  solely  as  a 
function  of  change  in  temperature.  Since  replacement  of 
crustal  material  is  suggested  for  the  Rio  Grande  Rift  and 
Basin  and  Range  provinces,  but  not  in  the  Colorado  Plateau, 
lateral  homogeneity  as  assumed  by  the  program  is  also 
fallacious . 

With  these  restrictions  in  mind,  the  program  can  still 
be  used  in  a  general  way.  The  heat  flow  observed  at  the 
surface  and  calculated  by  the  model  does  not  rely  on  density 
changes  from  temperature  change  or  convection.  Because  no 
assumptions  were  violated,  the  program  is  considered  to  be 
an  accurate  predictor  of  near-surface  temperature  gradients. 
Although  assumptions  behind  uplift  appear  to  be  invalid  in 
New  Mexico  and  Arizona,  regional  uplift  of  the  area  is  well 
documented  and  is  reviewed  in  an  earlier  section  of  this 
report.  The  thinning  of  the  crust  and  introduction  of 
basalt  at  its  base  should  result  in  relative  downwarp  with 
respect  to  neighboring  regions,  and  is  observed  between  the 
Colorado  Plateau  and  the  extensional  provinces.  Since  this 
is  spatially  associated  with  high  heat  flow,  the  uplifts 
predicted  by  the  paradigms  in  this  study  should  be 
considered  as  maximum  values. 
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It  is  therefore  concluded  from  heat  flow  and  uplift 
comparisons  that  southeastern  Arizona  and  southwestern  New 
Mexico  is  the  site  of  two  thermal  anomalies:  A  long 
standing  (about  28  Ma)  thermal  pulse  at  the  base  of  the 
lithosphere  along  with  a  shorter  term  (5.4  Ma)  anomaly  at 
the  base  of  the  crust.  These  perturbations  are  likely  to  be 
related,  and  can  explain  the  present  sta  of  volcanism, 
uplift,  extension  and  rifting  in  southeastern  Arizona  and 
southwestern  New  Mexico.  Because  the  uplift  is  regional,  it 
is  reasonable  to  expect  that  the  deep  anomaly  is  broader 
than  the  shallow  pulse.  Also,  because  of  the  narrow  band  of 
high  heat  flow  and  recent  volcanic  activity,  the  shallow 
anomaly  is  inferred  to  be  spatially  associated  with  the 
southeastern  margin  of  the  Colorado  Plateau. 


CONCLUSIONS 


This  study  addresses  two  fundamental  problems  in 
southeastern  Arizona  and  southwestern  New  Mexico:  the 
present  day  crustal  geometry  and  the  post-Laramide  tectonic 
history  of  the  southeastern  margin  of  the  Colorado  Plateau. 

A  solution  to  the  former  problem  comes  from  the  combined 
analysis  of  seismic  refraction  and  gravity  data,  and  the 
latter  problem  is  investigated  through  the  use  of  present 
day  structure,  surface  geologic  information  and  heat  flow. 

To  date,  available  data  indicate  a  complex  history  that 
includes  compression,  uplift,  tension,  thermal  disturbances 
and  extensive  volcanism  throughout  the  region. 

Keller  and  others  (1975)  determined  that  the 
physiographic  and  tectonic  boundaries  between  the  Colorado 
Plateau  and  Basin  and  Range  provinces  do  not  coincide  in 
Utah.  The  boundary  is  therefore  thought  to  be  transitional 
instead  of  abrupt.  This  distinction  is  also  supported  in 
Arizona  (Braumbaugh,  1987).  The  absence  of  a  physiographic 
boundary  east  of  the  Mogollon  Rim  is  therefore  irrelevant  to 
the  determination  of  the  tectonic  transition  zone  in  New 
Mexico.  A  complicating  factor  in  the  investigation  of  the 
transition  zone  in  southeastern  Arizona  and  southwestern  New 
Mexico  is  the  presence  of  the  mid-Tertiary  Mogollon-Datil 
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volcanic  field  (Elston  and  others,  1968,  1970,  1976a).  The 
volcanics  hide  previous  features  and  have  influenced 
structural  response  to  subsequent  tectonic  activity. 

Because  of  the  volcanic  cover,  it  is  possible  to 
determine  crustal  structure  only  through  the  analysis  of 
geophysical  data  constrained  by  geologic  information. 
Previous  geologic  work  has  shown  that  a  cluster  of  calderas 
that  vented  exclusively  high-silica  rhyolite  is  found  on  the 
Mogollon  plateau  (see  e.g.  Elston  and  others,  1976a) . 
Although  two  calderas  of  the  high-silica  rhyolite  suite 
occur  outside  this  cluster,  they  are  volumetrically  minor 
with  respect  to  the  plateau.  Rhodes  (1976)  has  determined 
that  rocks  from  each  cauldron  have  certain  unique  chemical 
characteristics,  but  that  the  entire  complex  is  intruded  by 
dikes  of  consistent  composition  that  have  similarities  to 
each  cauldron  fluid.  Rhodes  called  the  rocks  contained  in 
these  dikes  "framework  lavas",  and  suggested  a  common  source 
for  all  of  the  volcanics.  Because  no  other  volcanic  suite 
was  vented  on  the  plateau  during  the  development  of  these 
calderas,  Elston  and  others  (1976a)  proposed  that  a  near¬ 
surface  magma  chamber  existed  at  that  time.  That  chamber 
was  thought  to  remain  as  a  pluton  (Coney,  1976a;  Krohn, 

1976) ,  and  was  restricted  to  the  area  below  the  Mogollon 
plateau . 
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The  integrated  use  of  seismic  refraction  and  gravity 
information  resolves  crustal  structure  to  an  extent  not 
possible  by  either  alone.  Analysis  of  the  data  indicate  the 
following: 

1)  Gravity  and  seismic  refraction  data  provides 
evidence  for  a  low-density,  low-velocity  body  in 
the  upper  crust  beneath  southwestern  New  Mexico 
and  southeastern  Arizona.  The  proximity  of  mid- 
Tertiary  volcanism  with  this  body  supports  the 
interpretation  of  the  presence  of  a  pluton  of 
similar  age  by  Elston  and  others  (1976a). 

However,  the  geophysical  data  from  both  profiles 
examined  in  this  study  indicate  that  the  pluton 
extends  about  50  km  north  of  the  boundaries 
suggested  by  Elston  and  others  (1976a) ,  Krohn 
(1976),  and  McIntosh  (1989). 

2)  The  boundary  of  the  pluton  as  predicted  by 
gravity  and  seismic  refraction  analysis 
corresponds  closely  to  the  -20  mgal  contour  line 
of  the  125  km  low-pass  gravity  map  over  the  study 
area.  If  this  contour  line  follows  the  outline  of 
the  pluton,  then  the  complex  is  larger  than  that 
predicted  by  the  profiles  used  in  this  study.  In 


fact,  it  would  be  approximately  230  km  by  155  km, 
extending  northwest  into  Arizona  (Figure  57) . 

3)  The  Moho  markedly  deepens  approximately  90  km 
north  of  Silver  City,  from  approximately  32  to  37 
km  over  a  distance  of  about  35  km.  This  region  is 
located  about  midway  between  the  north  and  south 
ends  of  the  plutonic  complex,  and  is  north  of  the 
high-silica  rhyolite  cauldrons  of  the  Mogollon 
plateau. 

4)  The  Moho  flattens  at  38-40  km  depth  about  150 
km  north  of  Tyrone.  This  point  is  defined  in  this 
study  as  the  southern  edge  of  the  Colorado 
Plateau.  The  final  model  describing  current 
crustal  structure  (Figure  38)  is  similar  to  that 
described  in  central  Arizona  by  Warren  (1969) . 
Therefore,  the  present  day  boundary  of  the 
Colorado  Plateau  trends  west-southwest  from 
central  Arizona  under  the  northern  edge  of  the 
Plutonic  complex,  then  swings  northward  along  the 
edge  of  the  Rio  Grande  Rift  (Figure  57) . 

5)  Gravi!:.y  data  indicate  that  the  Moho  may  become 
shallower  by  about  2  km  near  Grants.  However,  the 
data  could  also  be  explained  by  expansion  of  upper 
crustal  layers  at  the  expense  of  the  lower  crust 
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in  that  region,  leaving  the  Moho  flat.  Seismic 
data  was  not  obtained  in  that  area,  yielding  lower 
resolution  of  major  crustal  boundaries. 

6}  If  crustal  thickness  is  used  as  a 
characteristic  feature  of  the  provinces  in  this 
region,  then  the  plutonic  complex  overlaps  the 
boundary  between  the  Colorado  Plateau  and  the 
extensional  provinces  to  the  south  and  east.  The 
plutonic  complex  is  spatially  associated  with  the 
present  day  transition  zone  in  Arizona  and  New 
Mexico. 

The  presence  of  a  pluton  with  greater  extent  than  the 
region  under  the  Mogollon  plateau  points  out  some 
ambiguities  in  previous  work.  For  example,  Elston  and 
others  (1976a)  suggest  that  the  outer  grabens  which  wrap 
around  the  Mogollon  plateau  are  caused  by  Basin  and  Range 
type  extension,  and  that  the  wrapping  is  caused  by 
resistance  of  the  pluton  to  through-going  faulting.  Since 
the  pluton  clearly  extends  to  the  north  of  the  Plains  of  St. 
Augustine,  this  interpretation  may  be  doubtful.  Rhodes 
(1976)  shows  that  the  outer  rim  of  the  Mogollon  plateau  is 
associated  with  framework  rhyolite,  probably  derived  from 
faults  tapping  the  magma  body  itself.  The  graben  faults  may 
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Figure  57.  Extent  of  the  proposed  upper  crustal  pluton  and 
its  relationship  to  the  125  km  low-pass  gravity 
field.  The  southeastern  margin  of  the  Colorado 
Plateau  is  shown  by  the  dashed  line. 
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be  associated  with  this  plumbing  system.  Furthermore,  the 
interior  faults  of  these  grabens  show  more  movement  than  the 
exterior  faults  (Elston  and  others,  1976a) .  The  spatial  and 
chemical  associations  of  these  graben/ fault  systems  may 
indicate  that  the  volcanic  complex  itself  is  the  causative 
factor,  as  explained  below. 

Coney  (1976a)  shows  the  Mogollon  plateau  to  consist  of 
a  central,  relatively  undeformed  region  that  has  subsided 
with  respect  to  its  margin.  Its  boundaries,  in  turn,  are 
uplifted  and  partly  composed  of  framework  rhyolite.  The 
"Basin  and  Range”  type  faulting  and  associated  grabens  lie 
outside  this  rim.  Thus,  a  possible  explanation  for  faulting 
on  the  perimeter  of  the  Mogollon  plateau  is  the  collapse  of 
the  entire  high-silica  rhyolite  volcanic  field  above  a  huge 
vented  magma  chamber.  The  uplifted  rims  and  associated 
faults  are  merely  a  giant  ring-fracture,  magma  filled 
system.  This  would  explain  fracturing  and  faulting  in  terms 
of  a  high-silica  rhyolite  portion  of  a  larger  plutonic 
complex,  as  opposed  to  faulting  associated  with  Basin  and 
Range  type  extension. 

The  Plains  of  St.  Augustine  may  be  of  different  origin 
or  have  a  different  history  than  the  other  rim  grabens. 
Chapin  (1971)  suggests  that  this  may  be  a  splay  of  the  Rio 
Grande  Rift.  Chapin  cites  the  presence  of  left-lateral 
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Oblique  slip  faults  at  its  northern  and  southern  boundaries, 
that  the  down- faulted  section  almost  completely  includes  the 
Datil  group,  and  that  structural  relief  is  at  least  1.1  km. 
Whether  this  is  a  branch  of  the  rift  is  not  determined  in 
this  study,  but  it  is  noteworthy  that  this  offset  is  minor 
compared  to  5000  m  of  movement  at  the  edge  of  the  rift  in 
the  nearby  Lucero  uplift  area  (Callender  and  Zilinski, 

1976) .  Also,  no  crustal  thinning  associated  with  this 
feature  is  indicated  from  the  present  study. 

The  spatial  relationship  between  the  high-silica 
rhyolite  pluton  and  the  rest  of  the  complex  cannot  be 
uniquely  determined  in  this  study.  Furthermore,  it  remains 
undetermined  whether  the  high-silica  fluids  come  from  a 
separate  source  area  or  are  chemically  related  to  the  rest 
of  the  complex.  Elston  and  others  (1976a,  1976b)  show 
tendencies  for  different  volcanic  suites  to  have  different 
average  chemical  characteristics  (such  as  Sr/  Sr  values 
and  types  of  mineralization) ,  but  overlaps  occur  in  these 
data  sets.  Separation  of  the  high-silica  rhyolite,  calc- 
alkalic  rhyolite  and  calc-alkalic  andesite  suites  into 
unique  source  reservoirs  is  a  step  that  remains  to  be 
proven. 

An  alternative  interpretation  to  a  unique  high-silica 
magma  chamber  is  one  in  which  the  high-silica  fluids  are 
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derived  from  a  larger  magma  chamber.  From  sedimentary 
evidence,  the  crust  under  the  Colorado  Plateau  in  New  Mexico 
was  low  relative  to  the  south  immediately  prior  to  and 
during  initial  stages  of  mid-Tertiary  volcanism  (Gather  and 
Johnson,  1986) .  The  structural  trend  has  subsequently 
reversed  (Pierce  and  others,  1979;  Morgan  and  Swanberg, 

1985) .  It  is  possible  that  the  most  felsic  components, 
including  volatiles,  worked  their  way  to  the  top  of  the 
chamber  as  it  developed.  These  pooled  under  the  Mogollon 
plateau  for  the  simple  reason  that  it  was  the  highest  part 
of  the  pluton  at  the  time.  Subsequent  cauldron-producing 
volcanic  events  were  above  the  lightest,  hottest  material. 

The  two  high-silica  cauldron  events  outside  the 
Mogollon  plateau  (Ratte  and  others,  1969;  Rhodes  and  Smith, 
1972)  may  reflect  local  maxima  in  the  roof  of  the  magma 
chamber.  The  density/cheroistry  stratification  would  satisfy 
the  requirement  that  only  high-silica  rhyolites  were  vented 
within  the  Mogollon  plateau  during  that  time.  This  pooling 
effect  may  have  caused  a  minor  local  uplift  over  the 
Mogollon  plateau,  creating  the  radial  tension  cracks 
reported  by  Elston  and  others  (1976a) .  Subsequent  collapse 
above  a  venting  magma  chamber  caused  normal  faulting  at  the 
perimeter  of  the  MP.  Rocks  of  the  calc-alkaline  suite  could 
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only  be  vented  from  the  south  end  of  the  pluton  after  the 
overlying  material  was  removed. 

The  cauldrons  which  did  not  vent  high-silica  rhyolite 
rocks  yet  are  associated  with  the  plutonic  complex  are  not 
as  clustered  as  those  on  the  Mogollon  plateau.  Rhyolites 
are  mapped  in  the  Datil  mountains  and  continue  north  to 
outcrops  of  Mesozoic  sediments  (New  Mexico  Geological 
Society  Map,  1982) .  Intermediate  volcanic  rocks  are  also 
observed  in  east  central  Arizona,  but  much  of  that  region  is 
covered  by  more  recent  basalt  flows  (Wilson  and  others, 

1969) .  However,  these  are  not  of  the  high-silica  type  found 
on  the  Mogollon  plateau  (Elston  and  others,  1976a) .  If 
these  rocks  were  also  vented  from  the  pluton,  a  chemical 
differentiation  may  have  existed  and/or  developed  during  the 
time  of  extrusion. 

Other  intermediate,  calc-alkaline  plutonic  events  in 
New  Mexico  and  Arizona  occur  outside  the  boundary  of  the 
complex  shown  in  Figure  57.  This  does  not  negate  the 
argument  that  the  magma  chamber  controlled  volcanism  where 
it  was  present.  The  entire  region  was  affected  by  calc- 
alkaline  volcanism  as  shown  in  Figure  10,  so  the  presence  of 
many  cauldrons  in  the  areas  is  expected  during  the  mid- 
Tertiary.  However,  this  interpretation  concludes  that  the 
high-silica  rhyolite  suite  developed  as  a  result  of  the 
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ability  to  tap  a  magma  chamber  with  a  maximum  volume  of 
approximately  180,000  km^. 

The  presence  of  a  late  Cenozoic  episode  of  crustal 
thinning  is  supported  in  the  Basin  and  Range  province  (Sinno 
and  others,  1981) .  Also,  low  Pn  velocities  in  the  Rio 
Grande  Rift  (Sinno  and  others,  1986)  and  in  this  study 
indicate  substantial  present-day  heating  at  the  base  of  the 
crust.  Recent  regional  uplift  commenced  around  5.5  Ma 
(Morgan  and  Swanberg,  1985) ,  and  includes  the  Colorado 
Plateau. 

Despite  this  heating  event,  the  crust  under  the 
Colorado  Plateau  appears  to  be  resistant  to  thinning  with 
respect  to  the  Basin  and  Range  province.  This  may  imply  a 
fundamental  chemical  difference  in  the  lower  crust  between 
the  Colorado  Plateau,  Rio  Grande  Rift  and  Basin-Range 
provinces.  If  present  day  crustal  geology  reflects  that  of 
mid-Tertiary  time,  then  perhaps  the  high-silica  rhyolite 
suite  is  a  product  of  this  difference  in  crustal  chemistry. 

A  predictive  model  that  calculates  heat  flow  and 
isostatic  uplift  has  been  constructed  to  determine  a 
temperature  history  of  the  region.  Within  the  limitations 
placed  on  this  model,  a  preliminary  conclusion  is  that 
present  day  heat  flow  and  elevation  of  the  surface  is  best 
explained  by  two  thermal  anomalies:  1)  a  thermal  pulse  that 
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requires  a  maximum  temperature  increase  of  400°C  above 
stable  continental  values  at  100  km  depth  has  existed  at 
least  since  mid-Tertiary  time.  This  has  resulted  in 
regional  isostatic  uplift  of  about  1.5  km.  2)  a  temperature 
perturbation  of  about  350°C  at  the  base  of  the  crust  was 
established  approximately  5.4  Ma  before  present. 

The  deep  anomaly  yields  a  1.5  km  uplift  at  steady- 
state.  Since  28  Ma  is  required  for  steady-state,  this 
implies  onset  of  the  thermal  pulse  during  the  mid-Tertiary 
as  a  minimum  estimate  if  temperatures  have  equilibrated. 
However,  volcanism  is  present  before  this  time  in  the 
region,  which  means  this  perturbation  may  be  older  than  the 
minimum  time  needed.  The  earliest  time  for  this  pulse  to  be 
initiated  is  marked  by  the  Cretaceous  Mancos  shale,  required 
to  be  at  or  below  sea  level  at  the  time  of  deposition 
(Morgan  and  Swanberg,  1985)  .  The  response  to  this  anomaly 
may  have  been  complicated  by  factors  such  as  lithospheric 
strength,  horizontal  and  vertical  stress. 

The  younger  event  is  likely  related  to  the  older, 
deeper  anomaly  and  corresponds  to  the  time  of  recent  uplift 
(Morgan  and  Swanberg,  1985) ,  rifting  (Seager  and  Morgan, 
1979)  and  volcanism  (Luedke  and  Smith,  1978) .  It  is 
spatially  associated  with  the  high  heat  flow  zone  at  the 
southeastern  boundary  of  the  Colorado  Plateau.  It  also 
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appears  to  be  consistent  with  present  day  crustal  structure 
determined  by  seismic  refraction  and  gravity  analysis  in 
this  study. 
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ABBREVIATIONS 


BRP 

CP 

Gsi 

GPP 

HFU 

Ma 

MDVF 

MP 

PCR 


Pg 

PgR 

PmP 

Pn 

RGR 

Sn 


Basin  and  Range  province 
Colorado  Plateau 
Giga-annum 

Great  Plains  province 
Heat  flow  unit  =  41.87  mWm’^ 

Mega -annum 

Mogollon-Datil  volcanic  field 
Mogollon  plateau 

P-wave  reflection  from  the  top  of  the  lower  crustal 
layer 

P-wave  refracted  through  the  upper  crust 

P-wave  reflection  from  the  top  of  the  intermediate 

crustal  layer 

P-wave  reflection  from  the  top  of  the  Moho 
P-wave  refracted  through  the  upper  mantle 
Rio  Grande  Rift 

S-wave  refracted  through  the  upper  mantle 
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LIST  OF  SOURCE  PARAMETERS 
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Notes: 

Base  station  coordinates  =  -3,439.369  3,611,737.316 

UTM  coordinates  referenced  to  -108.5°  longitude 

*  Origin  times  to  the  nearest  minute:  at  least  one  record 
from  these  dates  occupied  the  same  site  as  a  record 
with  a  positively  identified  origin  time. 
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LIST  OF  RECORDING  STATIONS,  LOCATIONS,  AND  ELEVATIONS 
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STATION 

TATITUDE 

T.ONGTTUDE 

QUADRANGLE  ELEVATION 

(FT) 

base 

32.659523 

108.341568 

Tyrone 

6425 

ts-1 

32.697453 

108.338661 

Tyrone 

5989 

ts-2 

32.712688 

108.296310 

Tyrone 

5700 

ts-3 

32.740116 

108.274460 

Tyrone 

5948 

tg-2 

32.841309 

108.223999 

Ft .  Bayard 

6730 

tg-3 

32.872486 

108.215538 

Ft .  Bayard 

6925 

tg-4 

32.894783 

108.234818 

Twin  Sister 

6570 

tg-5 

32.909691 

108.228325 

Twin  Sister 

6760 

tg-7 

32.945950 

108.196449 

Twin  Sister 

7475 

tg-9 

32.975708 

108.207985 

Twin  Sister 

7250 

tg-10 

33.015926 

108.229034 

Copperas  Peak 

6920 

tg-11 

33.039646 

108.219048 

Copperas  Peak 

5840 

tg-12 

33.061916 

108.209946 

Copperas  Peak 

6120 

tg-13 

33.098026 

108.194870 

Copperas  Peak 

6640 

tg-14 

33.120518 

108.189758 

Copperas  Peak 

7400 

tg-15 

33.159019 

108.197685 

Gila  Hot  Spr. 

6440 

tg-16 

33.202793 

108.214050 

Gila  Hot  Spr. 

5600 

tg-17 

33.185028 

108.028549 

Middle  Mesa 

6780 

tg-18 

33.214706 

108.028084 

Middle  Mesa 

7520 

tg-19 

33.241726 

108.061966 

Middle  Mesa 

7080 

tg-20 

33.275455 

108.060669 

Wall  Lake 

6510 

tg-21 

33.307198 

108.070786 

Wall  Lake 

6910 

tg-22 

33.340206 

108.063293 

Wall  Lake 

6510 

tg-2  3 

33.373352 

108.080505 

Wall  Lake 

7090 

tg-27a 

33.450386 

108.054703 

Spring  Canyon 

7000 

tg-28a 

33.478180 

108.017639 

Taylor  Peak 

7118 

tg-29a 

33.506470 

107.980049 

Indian  Peaks  E. 

7204 

tg-30a 

33.534870 

107.935738 

Indian  Peaks  E. 

7222 

tg-30 

33.547806 

107.919258- 

Indian  Peaks  E. 

7300 

tg-31 

33.566269 

107.905579 

Indian  Peaks  E. 

7325 

tg-32 

33.595684 

107.883057 

Indian  Peaks  E. 

7398 

tg-3  3 

33.631451 

107.872169 

Paddy's  Hole 

7460 

tg-34 

33.659649 

107.858421 

Paddy's  Hole 

7515 

tg-3  6 

33.722771 

107.820694 

Paddy's  Hole 

7600 

tg-3  8 

33.760544 

107.790253 

Luera  Mts .  E . 

7460 

tg-42a 

33.896481 

107.832504 

C-N  Lake 

7105 

tg-43a 

33.911167 

107.844986 

C-N  Lake 

7045 

tg-44a 

33.930710 

107.848915 

C-N  Lake 

7060 

tg-4 5a 

33.947365 

107.840919 

C-N  Lake 

6965 

tg-46a 

33.966145 

107.835365 

C-N  Lake 

7090 

tg-47a 

33.991222 

107.823097 

C-N  Lake 

7132 

tg-4 8a 

34.031307 

107.808899 

Anderson  Pk. 

7079 

tg-49a 

34.056866 

107.805618 

Anderson  Pk. 

7156 

tg-50a 

34.069225 

107.773674 

Anderson  Pk. 

7063 

tg-51a 

34.091316 

107.774948 

Anderson  Pk. 

7119 
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STATION 

LATITUDE 

LONGITUDE 

QUADRANGLE  ELEVATION 

(FT) 

tg-43 

33.908760 

107.781120 

C-N  Lake 

6910 

tg-45 

33.986760 

107.800285 

C-N  Lake 

6897 

tg-46 

34.010395 

107.774895 

Anderson  Pk. 

6968 

tg-47 

34 . 039173 

107.768524 

Anderson  Pk. 

6999 

tg-48 

34.067383 

107.768425 

Anderson  Pk. 

7050 

tg-49 

34.093628 

107.768356 

Anderson  Pk. 

7096 

tg-50 

34 . 169605 

107.808144 

Datil 

7519 

tg-52 

34.207268 

107.862831 

Datil 

7590 

tg-54 

34.239124 

107.831161 

Datil 

7800 

tg-55 

34.264240 

107.810883 

Cal  Ship  Mesa 

7940 

tg-56 

34.274857 

107.789024 

Cal  Ship  Mesa 

8120 

tg-58 

34.307701 

107.759628 

Cal  Ship  Mesa 

8100 

tg-59 

34.330269 

107.767693 

Cal  Ship  Mesa 

7780 

tg-58a 

34.305569 

107.688385 

Dog  Springs 

7361 

tg-59a 

34.327656 

107.713844 

Dog  Springs 

7590 

tg-60 

34.344658 

107.743675 

Dog  Springs 

7520 

tg-61 

34.359093 

107.725433 

Dog  Springs 

7490 

tg-62a 

34 . 383461 

107.683365 

D  Cross  Mtn. 

7365 

tg-63a 

34.413815 

107.667313 

D  Cross  Mtn. 

7180 

tg-64 

34 . 442356 

107.653877 

D  Cross  Mtn. 

6760 

tg-65 

34.468246 

107.637146 

D  Cross  Mtn. 

6482 

tg-66 

34.493443 

107.577156 

Table  Mountain 

6320 

tg-67 

34.515568 

107.578476 

Pueblo  Viejo  M. 

6430 

tg-68 

34.536659 

107.579109 

Pueblo  Viejo  M. 

6470 

tg-69 

34.554462 

107.596077 

Pueblo  Viejo  M. 

6540 

tg-70 

34.583324 

107.571724 

Pueblo  Viejo  M. 

6630 

tg-74 

34.689655 

107.596329 

Broom  Mountain 

7260 
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LIST  OF  RECORDING  STATIONS  AND  DISTANCE  TO  EACH  SITE 
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STATION 

DISTANCE 

ts-1 

4.40 

ts-2 

7.30 

ts-3 

11.22 

tg-1 

20.40 

tg-2 

22.97 

tg-3 

26.40 

tg-4 

27.91 

tg-5 

29.76 

tg-6 

31.88 

tg-7 

34.50 

tg-8 

35.54 

tg-9 

37.20 

tg-10 

40-92 

tg-11 

43.72 

tg-12 

46.30 

tg-13 

50.59 

tg-14 

53.07 

tg-15 

57.02 

tg-15.5 

59.07 

tg-16 

61.44 

tg-17 

65.18 

tg-18 

68.16 

tg-19 

69.63 

tg~20 

73.16 

tg-21 

76.14 

tg-22 

79.86 

tg-23 

82.83 

tg-24 

85.71 

tg-25 

89.90 

tg-27a 

91.71 

tg-2  6 

94.00 

tg-28a 

95.68 

tg-27 

97.49 

tg-28 

99.75 

tg-2 9a 

99.79 

tg-2  9 

102.47 

tg-30a 

104.21 

tg-30 

106.09 

tg-31 

108.48 

tg-32 

112.30 

tg-3  3 

116.34 

tg-34 

119.71 

tg-35 

122.71 

tg-3  6 

127.52 

STATION 

DISTANCE  fKM^ 

tg-37 

129.99 

tg-38 

132.46 

tg-39 

137.27 

tg-42a 

145.13 

tg-42 

145.5 

tg-43a 

146.31 

tg-43 

148 . 06 

tg-44a 

148 . 24 

tg-45a 

150.21 

tg-46a 

152.37 

tg-44 

152.50 

tg-44.5 

153.42 

tg-47 

155.36 

tg-45 

155.58 

tg-4 5. 5 

156.94 

tg-4  6 

158.84 

tg-48a 

159.99 

tg-46.5 

160.63 

tg-47 

162.01 

tg-49a 

162.78 

tg-48 

164.98 

tg-50a 

165.03 

tg-51a 

167.33 

tg-49 

167.79 

tg-49.75 

172.96 

tg-50 

174.63 

tg-50 . 5 

175.86 

tg-52 

177.32 

tg-53 

179.44 

tg-54 

181.48 

tg-55 

184.68 

tg-56 

186.35 

tg-57 

188.82 

tg-57a 

190.48 

tg-58 

190.62 

tg-58a 

192.37 

tg-59 

192.80 

tg-59a 

193.98 

tg-60 

194.95 

tg-61 

197.00 

tg-62a 

200.71 

tg-63a 

204.38 

tg-64 

207.77 

tg-65 

210.97 
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DISTANCE 


STATION  DISTANCE 


tg-66 

215.40 

tg-67 

217.68 

tg-68 

219.87 

tg-69 

221.25 

tg-70 

224.98 

tg-71  227.38 
tg-72  229.48 
tg-73  233.41 
tg-74  235.51 
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PROCEDURES  FOR  SEISMIC  DATA  RETRIEVAL 
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Flow  Chart  of  Seismic  Processing 


Cubic  spline  interpolation;  record  section 
plotting  for  subsequent  analysis 


APPENDIX  V 


PROGRAM 


HF2D 
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c********************************************************** 

c  * 

c  Program  HF2D  by  R.  V.  Schneider  * 


c  Program  to  determine  the  flow  of  heat  in  a  * 
c  two-dimensional  region  using  the  Alternating  * 
c  Direction  Implicit  (ADI)  method  * 


c 

* 

c 

This  model  assumes  that  the  number  of  nodes  in 

* 

c 

the  X- 

and  z-directions  are  equal.  Also,  the 

* 

c 

side  boundary  conditions  are  that  lateral  heat 

* 

c 

flow  is 

set  equal  to  zero,  i.e.  all  heat  flow 

★ 

c 

is  in  the  vertical  direction. 

* 

u 

c 

This  program  requires  the  following  data  for 

* 

c 

setup: 

* 

■k 

c 

Card  1; 

Format (is, 2 f 5. 1) 

* 

c 

Isize 

Number  of  nodes  on  each  side  (square 

* 

c 

region,  max  101)  (NOTE:  Isize 

* 

c 

must  be  an  ODD  NUMBER) 

•k 

c 

Delx 

Distance  (km)  between  nodes  in  the 

* 

c 

x-direction 

* 

c 

Delz 

Distance  (km)  between  nodes  in  the 

* 

c 

z-direction 

* 

c 

c 

Card  2: 

Format (3 is) 

* 

c 

Idt 

Time  increment  (thousands  of  years) 

* 

c 

between  time  steps 

* 

c 

Numt 

Maximum  number  of  time  steps  in  model 

* 

c 

evolution 

* 

c 

Itstep 

Values  for  all  parameters  will  be 

* 

c 

written  to  file  at  increments  of  Itstep. 

* 

c 

E.G.  for  Idt  of  10  and  Itstep  of  SO, 

* 

c 

parameter  values  will  be  written  every 

* 

c 

S00,000  years 

* 

k 

c 

c 

Card  3: 

Format (2f 6. 5) 

k 

c 

Diff 

Diffusivity  of  material  (cm  squared  per 

k 

c 

second) 

k 

c 

Cond 

Conductivity  ,  (W  /cm  C) 

k 

k 

c 

Card  4 : 

Format (2fS. 1, f6. 4) 

k 

c 

Tt 

Top  temperature  of  model 

k 

c 

Coef 

Leading  coefficient  for  a  gaussian 

k 

c 

temperature  distribution  at  the  bottom 

k 

c 

boundary  of  the  model 

k 

Coefficient  of  the  power  term  in  the 
gaussian  equation 
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c********************************************************** 

c 

Program  hf2d( infile) 
c 

Integer  Isize, Idt ,Numt, Itstep,Num,maxdim 

Parameter  (maxdim=101) 

Real  T(maxdim),  Temp(maxdim,maxdim) , 

&  Last (maxdim,maxdim) , 

&  Temp2 (maxdim,maxdim) ,ba(maxdim) , 

&  a (maxdiro) , b (maxdim) , c (roaxdim) , 

&  r (maxdim) ,u (maxdim) ,Q (maxdim) ,Rho (maxdim) , 

&  Tempi (maxdim, maxdim) iu (maxdim) 

Real  Alpha, Beta, Amp, Lamp, Coef, Pc, Cond, 

&  Delx, Delz, Delt,Tt 

Character*20  infile 

Common  /parameters/  Isize, Texp, Cond, Delx, Delz 
c 

0pen(10, file=infile,status=*old' ) 

Read (10, 6)  Isize, Delx, Delz 
Read (10, 3)  Idt,Numt, Itstep 
Read (10, 4)  Diff,Cond 
Read (10, 5)  Tt,Coef,Pc 
Close(lO) 

3  Format ( 3 i5) 

4  Format (2f 6. 5) 

5  Format(2f5. 1, f6.4) 

6  Format (i5,2f 5.1) 
c 

Delt  =  Float (Idt)  *  1000. 

Isize  =  (Isize/2)  *2+1 
D  =  Diff  *  .003156 
Cond  =  Cond  *  l.e5 

Alpha  =  (D  *  Delt)  /  (  2.  *  (Delx)**2  ) 

Beta  =  (D  *  Delt)  /  (  2.  *  (Delz)**2  ) 

Texp  =  .00003 
Num  =  0 
c 

c********************************************************** 

c  * 

c  Set  up  the  initial  model  * 

c  * 

c********************************************************** 

c 

c*****  Initial  temperature  distribution 
c 

Call  tempinit (Isize , delz , tt , T , maxdim) 
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Do  8  i=l,Isize 
Q(i)  =  0.0 
ba ( i )  =  0.0 
iu(i)  =  0.0 
Do  7  j=l,Isi2e 
Teinp(i,j)  =  T(j) 

Teinp2(i,j)  =  T(j) 

7  Continue 

8  Continue 
c 

c*****  Temperature  boundary  condition,  bottom 
c 

Do  15  i=l,Isize 

Temp(i, 1)=T(1) +Coef*exp(- (Pc* (Float (i-1) *Delx) **2) ) 
Temp2 (i, 1) =T ( 1) +Coef*exp (- (Pc* (Float (i-1) *Delx) **2) ) 
Temp(i,Isize)  =  Tt 
Temp2 (i,Isize)  »  Tt 

15  Continue 
c 

c*****  Initial  density  distribution 
c 

Call  density ( Isize , delz , Rho , T , maxdim) 
c 

c*****  Write  initial  values  to  file 
c 

Open ( 11 , f ile= ' rhoO ' , status= ' new ' ) 

Open ( 13 , f ile= ' baO • , status= ' new ' ) 

Open ( 14 , f ile= ' iuO ' , status^ ' new ' ) 

Open ( 15 , f ile= ' TO ' , status= ' new ' ) 

Do  11  j=l, Isize 

Write (11, 13)  j,rho(j) 

Write(13,13)  j,ba(j) 

Write (14, 13)  j,iu(j) 

Write(15,13)  j ,T(Isize-j+l) 

11  Continue 

13  Format ( i3, 2x, f8. 3) 

14  Format (2i3 , 2x, f7 . 2) 

20  Format (11 (lx, f 6, 1) ) 

Close (11) 

Close (13) 

Close (14) 

Close (15) 

Open ( 12 , f ile= ' tmO • , status* ' new ' ) 

Do  22  j=l, Isize 
Do  21  i=l, Isize 

Write(12,14)  i, j ,Temp(i, j) 

21  Continue 

22  Continue 
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Close (12) 
c 


c  * 
c  Begin  the  calls  to  the  subroutine  'tridag'  for  * 
c  heat  flow  evolution.  All  of  this  will  be  inside  * 
c  of  a  large  loop  that  will  increment  time  over  * 
c  values  set  by  the  user  * 
c  * 


Ck -k  ic  icic -k  it -k -k  ic  *  ic -k -k -k -k -k  It  it -k  It -k  *  ic  ** -k  it  if  It  4e  ** -k  it -k  ie  **  it  * -k  ic  it -k  * -k -k -k -k -k  is  *  ic  **  1:  it 

C 

Do  999  k  =  l,Numt 
c 

b(l)  =  1.  +  (2.  *  Alpha) 
c(l)  =  -2.  *  Alpha 
b(Isi2e)  =  b(l) 
a(Isize)  =  c(l) 
c 

Do  70  n=2,Isize-l 
a(n)  =  -Alpha 
b(n)  =1.  +  (2.  *  Alpha) 
c(n)  =  -Alpha 
70  Continue 
c 

c*****  The  first  call  to  tridag  involves  holding  j 
c*****  constant  while  incrementing  i.  That  is,  each  row 
c*****  will  be  calculated  for  the  first  call,  the  '1/2 
c*****  Time  step', 
c 

Do  200  j  =  2,Isize-l 
c 

c*****  Determine  RHS  values  before  calling  tridag 
c 

Do  100  i=l,Isize 

r(i)  =  Beta  *  (Temp(i,j-1)  +  Temp(i,j+1))  + 

&  (1.  -  (2.  *  Beta))  *  Temp(i,j) 

100  Continue 

c 

c*****  First  call  to  tridag 
c 

Call  tridag(a,b,c,r,u, Isize,maxdim) 
c 

c*****  Now  place  the  u  vector  in  its  correct  spot  in 
c*****  the  new  model: 
c 

Do  150  i=l,Isize 

Temp2(i,j)  =  u(i) 

150  Continue 
200  Continue 
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C 

204 

c 

c***** 

c 


250 

c 

c***** 

c***** 

c***** 

c 


& 

270 


& 

290 


& 

310 


C 


350 

400 

c 

c 


Format ( 11 ( lx, f 6.1) ) 

Set  up  a,b,c's  for  the  second  call 

b(l)  =  1. 
c(l)  =  0. 
a(Isize)  =  0. 
b(Isize)  =  1. 
r(Isize)  =  Tt 
Do  250  m=2,Isize-l 
a(m)  =  -Beta 
b(m)  =  1.  +  (2.  *  Beta) 
c(m)  =  -Beta 
Continue 

Now  set  up  the  second  call  by  holding  i  constant 
and  incrementing  j .  Thus  the  columns  will  now 
be  calculated. 

Do  400  i=l,Isize 
r(l)  =  Temp2(i,l) 

If(i.eq.l)  Then 

Do  270  j*2,Isize-l 

r(j)  =  (1.  -(2.*Alpha))*Temp2(l,j)  + 

(2 . *Alpha) *Temp2 (2 , j ) 

Continue 

Else  If (i.eq.Isize)  Then 
Do  290  j=2,Isize-l 

i’(j)  =  (!•  -(2.*Alpha)  )*Temp2(Isize,  j) 

+  (2 . *Alpha) *Temp2 (Isize-1, j ) 

Continue 

Else 

Do  310  j=2,Isize-l 

r(j)  -  Alpha*(Temp2 (i+i, j)+Temp2 (i'l, j) )  + 
(1.-  2 .*Alpha) *Temp2 (i, j ) 

Continue 

Endif 

Call  tridag (a,b,c,r,u,Isize, maxdim) 

Do  350  j=2, Isize-1 
Temp(i,j)  =  u(j) 

Continue 

Continue 

Check  evolution  of  model  after  this  time  step 


Lamp  «  0.0 
Do  600  i«l,Isize 
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Do  550  j=l,Isize 

Amp  =  Abs(Temp(i, j)  -  Last(i,j)) 

If (Amp. gt. Lamp)  Lamp  =  Amp 
550  Continue 
600  Continue 

If (Lamp. le. 0.1)  go  to  1000 
Do  800  i=l,Isize 
Do  750  j=l,Isize 

Last(i,j)  =  Temp(i,j) 

750  Continue 
800  Continue 
c 

c*****  Determine  if  output  file  should  be  written 
c 

If  (mod(lc,Itstep)  .eq.O)  Call 
&  evolve (T, Temp, Rho,Num,maxdim) 

c 

999  Continue 

1000  Call  evolve (T, Temp, Rho,Num,maxdim) 
l=k-l 

Open ( 15 , f ile= • tstep ' ) 

Write (15,*)  'Number  of  time  steps  =  ',1 
Close(15) 

Stop 

End 

c 

Subroutine  tempinit ( Isize , delz , tt , T , maxdim) 
c 

c*********************************************************** 
c  * 

c  Subroutine  to  set  up  a  typical  temperature  * 

c  distribution  within  the  lithosphere  * 

c  * 

c*********************************************************** 
c 

Real  qs , qm , tt , hg , tcon , depth , del z , Ten 
Dimension  T (maxdim) 

Integer  Isize 
c 

c*****  Set  up  initial  parameters 
c 

qs  =  15000. 
qm  =  7000. 
hg  *  800. 
tcon  *  600. 
c 

Do  10  j=l, Isize 

depth  =  float (j-1)  *  delz 
If  (depth. le. 10.0)  Then 
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T(Isize-j+l)  =  tt  +  (qs  *  depth) /tcon- (hg  * 

&  depth**2)/ (2 . *tcon) 

Ten  =  T(Isize-j+i) 

Else 

T(Isize-j+i)  =  Ten  +  qm  *  (depth  -  10.)  /  tcon 
End  if 

10  Continue 
Return 
End 
c 

Subroutine  density (Isize, delz ,Rho,T,inaxdiin) 
c 

c*********************************************************** 


c  * 
c  Subroutine  to  initialize  density  values  * 
c  in  the  lithosphere  * 
c  * 


Ck-khie-k -k-k  *1i-k -k -kit  1e*-k-k-k -k  * -kick  * -kit*  *11******  *****  ** -kick  Ic*  * -kit -k  * -k  *  ic  * -k  ** 
C 

Real  Texp,delz 

Dimension  Rho (maxdim) ,T(maxdim) 

Integer  Isize 
Texp  =  3.E-5 
,c 

Do  25  j=l, Isize 

Depth  =  Float (j-1) *delz 
If (Depth. le. 40.0)  Then 

Rho(j)  =  2,6  +  (0.01*depth) 

Else 

Rho(j)  =  2.3  +  (1.  -  Texp  *  (T(j)  -  T(l))) 

Endif 

25  Continue 
Return 
End 
c 

Subroutine  tridag (a,b,c,r,u,n, maxdim) 
c 

Ck**************************************************** ****** 


C  * 

c  Solves  for  a  vector  u  of  length  n,  the  tridagonal  * 
c  linear  set  given  by  equation  (2.6.1),  in  * 
c  "Numerical  Recipes  in  Fortran" .  * 
c  * 
c  a,  b,  c,  r  are  input  vectors  and  are  not  modified.  * 
c  * 


C*  *  "k  *********  k  ***********************  k*  ***************  ****** 
C 

Dimension  gam (maxdim) , a (maxdim) ,b (maxdim) ,c (maxdim) , 

&  r (maxdim) ,u (maxdim) 
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Real  bet 
c 

If  (b(l).eq.O.)  Pause 
bet  =  b(l) 
u(l)  =  r(l)/bet 
Do  50  j=2,n 

gain(j)  =  c(j-l)/bet 

bet  =  b(j)  -  a(j)  *  gain(j) 

If (bet. eq. 0 . )  Pause 
u(j)  =  (r(j)  -  a(j)  *  u(j-l))/bet 
50  Continue 

Do  60  j=n-l,l,-l 

u(j)  =  u(j)  -  gain(j+l)  *  u(j+l) 

60  Continue 
Return 
End 
c 

Subroutine  evolve (T , Temp , Rho , Hum , maxdim) 
c 

c********************************************************** 

c  * 

c  Subroutine  to  calculate  output  values  * 

c  * 

c***********************^********************************** 
c 

Integer  maxdim 

Dimension  T (maxdim) , Temp (maxdim, maxdim) , 

&  Rho (maxdim) ,F (maxdim) 

Real  hf (maxdim) ,ba (maxdim) , Drho (maxdim, maxdim) , 

&  iu (maxdim)  Sumout,out, 

&  DDrho (maxdim*2 , maxdim) ,TD (maxdim, maxdim) 

Integer  Isize,Num 
Character  str*4,Name*6 

Common  /parameters/  Isize,Texp, Cond, Delx, Delz 
c 

Num  =  Num  +  1 
Write (str, ' (i4) ' )  Num 
Call  STRIP (str) 
c 

c*****  Surface  heat  flow; 
c 

Name  =  'hfV/str 

Open ( 1 1 , f i 1 e=Name , status= ' new ' ) 

Do  10  i=l,Isize 

hf(i)  =  -Cond  *  (Temp(i,Isize)  - 
&  Temp(i,Isize-l) )  /  Delz 

hf(i)  =  hf(i)  /  4.l87e4 
Write(ll,lll)  hf(i) 
c  write (11, 113)  i,hf(i) 
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10  Continue 
Close(ll) 
c 

c*****  Establish  delta  rho 
c 

Do  15  i=l,Isize 
Do  14  j=l,Isize 

TD(i,j)  =  Teinp(i,j)  -  T(j) 

If ( j . eq. Isize)  Then 
Drho ( i , j )  =  0.0 
Else 

Drho(i,j)  =  -Rho(j) *Texp*TD(i, j) 

End  if 

14  Continue 

15  Continue 
c 

c*****  Isostatic  uplift: 
c 

Name  =  'iu'//str 

Open ( 12 , f ile=Name , status= ' new ' ) 

Do  19  i  =  1, Isize 
Do  18  j  =  1, Isize 
F(j)  “  Drho(i,j) 

18  Continue 

Call  Simpson(F,out,maxdim) 

iu(i)  =  -out 

Write (12, 111)  iu(i) 

19  Continue 
Close(12) 

c 

c*****  Write  temperature  distribution  to  file 
c 

Name  =  'tm'/Zstr 

Open ( 14 , f ile=Name , status* ' new ’ ) 

Do  65  j=l, Isize 
Do  55  i*l, Isize 

Write(14,114)  i, j,Temp(i, j) 

55  Continue 
65  Continue 
Close (14) 
c 

111  Format (f 15. 10) 

112  Format (i3,2x,f 15. 10) 

114  Format(2(2x,f4.0) ,2x,f7.2) 

115  Format(ll(lx,f6.1) ) 

116  Format(6(fl0.5,2x) ) 
c 

Return 

End 
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Subroutine  Siinpson(F,Suin,inaxdim) 
c 

c* ********************************************************** 

c  * 

c  Subroutine  to  apply  the  extended  Simpson's  rule  * 

c  procedure  * 

c  * 

c*********************************************************** 

c 

Real  Sum,Suinl,Sum2/delz 
Dimension  F(maxdim) 

Integer  Isize 

Common  /parameters/  Isize, Texp, Cond, Delx, Delz 
c 

Suml  =  (17.  *  (f(l)  +  f (Isize)) 

&  +  59.  *  (f(2)  +  f(Isize-l)) 

&  +  43.  *  (f(3)  +  f(Isize-2)) 

&  +  49.  *  (f(4)  +  f (Isize-3) ) )/48. 

Sum2  =  0.0 
Do  40  j=5,Isize-4 
Sum2  =  Sum2  +  f ( j ) 

40  Continue 

Sum  =  (Suml  +  Sum2)*delz 

Return 

End 
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